You need to enable JavaScript to run this app.
优惠活动
大模型
产品
解决方案
定价
更多

将Fortran/Julia的Sn尺度统计量O(nlogn)算法转Python遇索引问题求助

实现Sn尺度统计量的O(nlogn)算法问题

我正尝试实现Christophe Croux与Peter J. Rousseeuw在《Time-Efficient Algorithms for Two Highly Robust Estimators of Scale》中描述的、用于计算Sn(一种尺度统计量)的O(nlogn)算法。该论文提供了Fortran代码,我还找到了一份Julia实现(function scaleS!)并尝试将其转换为Python代码,但代码无法运行。我认为问题出在索引混淆上——Fortran和Julia的数组索引从1开始,而Python从0开始。以下是我的Python代码:

def sn(x):
    N = len(x)

    x.sort()
    a2 = [0] * N
    a2[0] = x[round(N / 2)] - x[0]

    for i in range(1, round((N + 1) / 2) - 1):
        nA = i - 1
        nB = N - i
        diff = nB - nA
        leftA = leftB = 1
        rightA = rightB = nB
        Amin = round(diff / 2) + 1
        Amax = round(diff / 2) + nA
        while leftA < rightA:
            length = rightA - leftA + 1
            even = 1 - (length % 2)
            half = round((length - 1) / 2)
            tryA = leftA + half
            tryB = leftB + half
            if tryA < Amin:
                rightB = tryB
                leftA = tryA + even
            elif tryA > Amax:
                rightA = tryA
                leftB = tryB + even
            else:
                medA = x[i] - x[i - tryA + Amin - 1]
                medB = x[tryB + i] - x[i]
                if medA >= medB:
                    rightA = tryA
                    leftB = tryB + even
                else:
                    rightB = tryB
                    leftA = tryA + even
        if leftA > Amax:
            a2[i] = x[leftB + i] - x[i]
        else:
            medA = x[i] - x[i - leftA + Amin - 1]
            medB = x[leftB + i] - x[i]
            a2[i] = min(medA, medB)

    for i in range(round((N + 1) / 2), N - 2):
        nA = N - i
        nB = i - 1
        diff = nB - nA
        leftA = leftB = 1
        rightA = rightB = nB
        Amin = round(diff / 2) + 1
        Amax = round(diff / 2) + nA
        while leftA < rightA:
            length = rightA - leftA + 1
            even = 1 - (length % 2)
            half = round((length - 1) / 2)
            tryA = leftA + half
            tryB = leftB + half
            if tryA < Amin:
                rightB = tryB
                leftA = tryA + even
            elif tryA > Amax:
                rightA = tryA
                leftB = tryB + even
            else:
                medA = x[i + tryA - Amin + 1] - x[i]
                medB = x[i] - x[i - tryB]
                if medA >= medB:
                    rightA = tryA
                    leftB = tryB + even
                else:
                    rightB = tryB
                    leftA = tryA + even
        if leftA > Amax:
            a2[i] = x[i] - x[i - leftB]
        else:
            medA = x[i + leftA - Amin + 1] - x[i]
            medB = x[i] - x[i - leftB]
            a2[i] = min(medA, medB)

    a2[N] = x[N] - x[round((N + 1) / 2)]
    lomed = a2[round((len(a2) + 1) / 2)]

    return lomed

内容的提问来源于stack exchange,提问作者Expertium

相关产品推荐
方舟 Agent Plan

超全模态模型 × Harness 升级,最新支持 Deepseek-V4.1-Flash、GLM-5.3 系列、Doubao-Seedream-5.0-pro、Kimi-K3 (部分), 限时 9.9 元起

最近更新时间:2026.08.01 08:56:10