我在努力实现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
1条答案
按热度按时间3mpgtkmj1#
下面是从Julia function直接翻译过来的代码,我通常使用的将Julia代码转换成Python的最快方法是:注解掉所有
end
,并将所有条件和循环放入:
,从所有索引中减去1
,最后替换等效函数。这些是代码中的一些特殊点,注意我删除了
lomed
函数,并添加了sum
,因为我不熟悉统计学中的低中值概念。0
开始索引数组,所以你需要调整代码中的一些索引来匹配它。div
,而Python的等价物是//
,所以需要将round(N / 2)
行更改为N // 2
。i
从1
到round((N + 1) / 2) - 1
,您需要调整切片x[i - tryA + Amin - 1]
中的索引,以便它们从0
而不是1
开始。i
从round((N + 1) / 2)
到N - 2
的第二个for循环中,您需要调整切片x[i + tryA - Amin + 1]
中的索引,以便它们从0
开始,而不是从1
开始。