基于Fortran/Julia代码的统计估计量的Python实现

umuewwlo  于 2023-02-11  发布在  Python
关注(0)|答案(1)|浏览(132)

我在努力实现O计算Sn的(nlogn)算法(规模的统计估计量)如Christophe Croux和Peter J. Rousseeuw的“Time-Efficient Algorithms for Two Highly Robust Estimators of Scale”中所述,本文提供了Fortran代码,我还找到了this Julia实现(function scaleS!),并试图将其转换为Python,但我的代码不起作用。我想这是因为我一直在搞乱索引-Fortran和Julia都是从1开始索引数组,Python从0开始索引。

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
3mpgtkmj

3mpgtkmj1#

下面是从Julia function直接翻译过来的代码,我通常使用的将Julia代码转换成Python的最快方法是:注解掉所有end,并将所有条件和循环放入:,从所有索引中减去1,最后替换等效函数。
这些是代码中的一些特殊点,注意我删除了lomed函数,并添加了sum,因为我不熟悉统计学中的低中值概念。

  • 因为Python从0开始索引数组,所以你需要调整代码中的一些索引来匹配它。
  • 因为Julia使用div,而Python的等价物是//,所以需要将round(N / 2)行更改为N // 2
  • 在for循环中,i1round((N + 1) / 2) - 1,您需要调整切片x[i - tryA + Amin - 1]中的索引,以便它们从0而不是1开始。
  • iround((N + 1) / 2)N - 2的第二个for循环中,您需要调整切片x[i + tryA - Amin + 1]中的索引,以便它们从0开始,而不是从1开始。
def sn(x):
    N = len(x)

    x = sorted(x)
    a2 = [0 for _ in range(N)]
    a2[0] = x[N//2]-x[0]
    for i in range(2,(N+1)//2+1):
        nA = i-1
        nB = N-i
        diff = nB-nA
        leftA = leftB = 1
        rightA = rightB = nB
        Amin = diff//2+1
        Amax = diff//2+nA

        while leftA < rightA:
            lent = rightA-leftA+1
            even = 1 - lent%2
            half = (lent-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-1]-x[i-tryA+Amin-2]
                medB = x[tryB+i-1] - x[i-1]
                if medA >= medB:
                    rightA = tryA
                    leftB = tryB+even
                else:
                    rightB = tryB
                    leftA = tryA+even
                # end
            # end
        # end
        if leftA > Amax:
            a2[i-1] = x[leftB+i-1]-x[i-1]
        else:
            medA = x[i-1]-x[i-leftA+Amin-2]
            medB = x[leftB+i-1]-x[i-1]
            a2[i-1] = min(medA,medB)
        #end
    #end

    for i in range((N+1)//2+1, N):
        nA = N-i
        nB = i-1
        diff = nB-nA
        leftA = leftB = 1
        rightA = rightB = nB
        Amin = diff//2+1
        Amax = diff//2+nA
        while leftA < rightA:
            lent = rightA-leftA+1
            even = 1-lent%2
            half = (lent-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]-x[i-1]
                medB = x[i-1]-x[i-tryB-1]
                if medA >= medB:
                    rightA = tryA
                    leftB = tryB+even
                else:
                    rightB = tryB
                    leftA = tryA+even
                # end
            # end
        # end
        if leftA > Amax:
            a2[i-1] = x[i-1]-x[i-leftB-1]
        else:
            medA = x[i+leftA-Amin]-x[i-1]
            medB = x[i-1]-x[i-leftB-1]
            a2[i-1] = min(medA,medB)
        # end
    # end # 20

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

    # Normalization
    cn = 1.0
    if N < 10:
        cn = [0, .743, 1.851, .954, 1.351, .993, 1.198, 1.005, 1.131][N-1]
    elif N%2 == 1:
        cn = N/(N-0.9)
    # end
    return cn * 1.1926 * sum(a2)

相关问题