2

我是一个非常新的 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)

所以 - 所有你在那里嗖嗖,你怎么想?谢谢!

4

1 回答 1

0

我查看了关于计算 Theil-Sen 斜率主题的R 文档,它们使用与 SciPy 相同的方法。

Conover (1980, p. 267) 建议使用以下截距估计器: 在此处输入图像描述

所以我猜 SciPy 方法很好。

于 2018-06-18T12:57:10.433 回答