将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
相关产品推荐
相关产品推荐

