Потенциальная ошибка Scipy в scipy.stats.mstats.theilslopes?

Question

Я довольно новый python / scipy / numpy и начал использовать его из-за встроенной функции оценки Theil-Sen в Scipy и дружественной итеративности Python. После сравнения результатов моего скрипта на python с другими вычислениями Theil-Sen, я обнаружил две ошибки в функции scipy.stats.mstats.theilslopes. Я надеюсь, что более опытный программист / статистик сможет подтвердить мои выводы.

Источник mstats ( https://github.com/scipy/scipy/blob/v0.14.0/scipy/stats/mstats_basic.py) имеет (я думаю) два раздела с ошибками. В первом разделе обе серии должны быть сделаны поплавками, и нет причин маскировать часть серии. Поэтому я бы пересмотрел этот код с:

  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()

Во-вторых, в вычислении перехвата Тейл-Сена, по-видимому, есть фундаментальная ошибка (как определено здесь: http://books.google.com/books?id=lK9gHXwYnqgC&pg=PA67). Текущий код вычисляет медиану для всех х и у, а затем работает на основе этих значений и наклона, чтобы получить перехват. Увидеть:

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

Я проверил R документацию на предмет расчета склонов Тейл-Сен, и там они используют то же самое, что и SciPy.

Коновер (1980, стр. 267) предлагает следующую оценку для перехвата:

Так что я думаю, что метод SciPy в порядке.