Scipy.stats.mstats.theilslopes中的潜在错误?

Question

我刚开始学习python、scipy和numpy，之所以使用它们是因为Scipy里有个内置的Theil-Sen估计器函数，而且Python的可迭代性很友好。在把我写的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()

改成：

  y = ma.asarray(y,dtype=float).flatten()

  n = len(y)
  if x is None:
      x = ma.arange(len(y), dtype=float)
  else:
      x = ma.asarray(x,dtype=float).flatten()

其次，在Theil-Sen截距的计算上似乎有一个根本性的错误（在这里定义：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)

那么——各位高手，你们怎么看？谢谢！

Answer 1

我查看了关于计算Theil-Sen斜率的R文档，发现他们使用的方法和SciPy是一样的。

Conover（1980年，第267页）建议使用以下公式来估算截距：

所以我觉得SciPy的方法是可以的。