我是一个非常新的python/scipy/numpy,并开始使用它是因为scipy内置的theilsen估计器函数和python友好的iterability。在将python脚本的结果与其他的Theil-Sen计算结果进行比较后,我认为我在scipy.stats.mstats.theilslopes函数。我希望一个更有经验的程序员/统计学家能证实我的发现。在
mstats源代码(https://github.com/scipy/scipy/blob/v0.14.0/scipy/stats/mstats_basic.py#L673)有两个部分(我认为)有错误。在第一部分中,两个系列都必须是浮动的,没有理由掩盖部分序列。因此,我将从以下方面修订本规范:
y = ma.asarray(y).flatten()
y[-1] = masked
n = len(y)
if x is None:
x = ma.arange(len(y), dtype=float)
else:
x = ma.asarray(x).flatten()
…收件人:
^{pr2}$第二,在计算泰尔森截距(定义如下:http://books.google.com/books?id=lK9gHXwYnqgC&pg=PA67#v=onepage&q&f=false)时,似乎有一个基本错误。当前代码计算所有x和y的中值,然后根据这些值和斜率计算截距。参见:
slopes = ma.hstack([(y[i+1:]-y[i])/(x[i+1:]-x[i]) for i in range(n-1)])
slopes.sort()
medslope = ma.median(slopes)
medinter = ma.median(y) - medslope*ma.median(x)
然而,正确的方法是将斜率应用于每个坐标对,然后根据这些值计算中值。所以,我认为正确的代码应该是:
slopes = ma.hstack([(y[i+1:]-y[i])/(x[i+1:]-x[i]) for i in range(n-1)])
slopes.sort()
medslope = ma.median(slopes)
intercepts = ma.hstack([(y[i] - medslope*x[i]) for i in range(n)])
intercepts.sort()
medinter = ma.median(intercepts)
所以——你们这些白痴,你们怎么看?谢谢!在
我检查了关于计算Theil-Sen斜率的R documentation,他们用的是SciPy。在
所以我想SciPy方法是好的。在
相关问题 更多 >
编程相关推荐