私はかなり新しい python/scipy/numpy であり、Scipy の組み込みの Theil-Sen 推定関数と Python のフレンドリーな反復可能性のために、それを使い始めました。私の python スクリプトの結果を他の Theil-Sen 計算と比較した後、scipy.stats.mstats.theilslopes 関数に 2 つの間違いを見つけたと思います。より経験豊富なプログラマー/統計学者が私の発見を裏付けることができることを願っています.
mstats ソース ( https://github.com/scipy/scipy/blob/v0.14.0/scipy/stats/mstats_basic.py#L673 ) には (私が思うに) 間違いのある 2 つのセクションがあります。最初のセクションでは、両方のシリーズを float にする必要があり、シリーズの一部をマスクする理由はありません。したがって、このコードを次のように修正します。
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)
それで - あなたがそこに飛び散るすべて、あなたはどう思いますか? ありがとう！
