2013-05-23 1 views
5

ブートストラップscikitsと回帰パラメータの信頼区間を推定私はちょうどscikitsを通じて利用できる素敵なブートストラップパッケージを試して始めました: https://github.com/cgevans/scikits-bootstrap(パイソン)

が、ために信頼区間を推定しようとしたとき、私は問題が発生しました線形回帰からの相関係数。返される信頼区間は、元の統計値の範囲外にあります。ここ

コードである:

import numpy as np 
from scipy import stats 
import bootstrap as boot 

np.random.seed(0) 
x  = np.arange(10) 
y  = 10 + 1.5*x + 2*np.random.randn(10) 
r0 = stats.linregress(x, y)[2] 

def my_function(y): 
    return stats.linregress(x, y)[2] 

ci = boot.ci(y, statfunction=my_function, alpha=0.05, n_samples=1000, method='pi') 

これは、CI = -0.605、0.644]の結果が得られるが、元の統計値は、R0 = 0.894です。

私はこれをRで試してみましたが、うまく動作しているようです。期待通りにciに跨っています。

助けてください!

答えて

8

Rコードを入力してください。

ここで問題となるのは、yをboot.ciに渡すだけですが、my_functionを実行するたびに、元のx( my_functionにx入力がないことに注意してください)。ブートストラッピングは再サンプリングされたデータに統計関数を適用するので、元のxyのサンプルを使用して統計関数を適用すると、無意味な結果になります。このため、BCAメソッドはまったく機能しません。実際には、同じ数の要素を持たないジャックナイフサンプルに統計関数を適用することはできません。

複数の1D配列を統計関数に渡す場合は、xy = vstack((x,y)).Tという複数の列を使用して、それらの列からデータを取り込むstat関数を使用できます。

def my_function(xysample): 
    return stats.linregress(xysample[:,0], xysample[:,1])[2] 

あなたはすべてのあなたのデータをいじって避けたかった場合あるいは、あなたはインデックスを操作する関数を定義することができ、その後、ちょうどboot.ciするインデックスを渡します

def my_function2(i): 
    return stats.linregress(x[i], y[i])[2] 

boot.ci(np.arange(len(x)), statfunction=my_function2, alpha=0.05, n_samples=1000, method='pi') 

。なお、いずれの場合も、BCAあなたが実際にパーセンテージ間隔を使いたい場合を除いて、method = 'bca'を使うこともできます。 BCAはずっと良いです。

私は、これらの方法の両方が理想的ではないことを認識しています。正直なところ、私はこのような複数の配列を自分のstat関数に渡す必要はありませんでした。そして、大部分の人々は、おそらくmeanをstatfunctionとして使用しています。私はここで最善のアイデアは、同じサイズの[0]配列のリストを渡すこと、例えばboot.ci([x,y],...)を許可し、それらのすべてを同時にサンプリングして、すべてをstatfunctionに別々の引数として渡すことです。その場合は、my_function(x,y)とすることができます。私はこれを行うことができるかどうかを見ていきますが、あなたのRコードを表示することができれば、それは素晴らしいことです。これに対処するより良い方法があるかどうかを見たいと思います。


更新:scikits.bootstrap(v0.3.1)の最新版において

、アレイの組を提供することができ、それらからサンプルをstatfunctionに別個の引数として渡されます。さらに、stat関数は配列出力を提供し、信頼区間は出力の各点について計算されます。したがって、これは現在非常に簡単です。以下はlinregressのすべての出力の信頼区間を与える:この場合

cis = boot.ci((x,y), statfunction=stats.linregress) 

cis[:,2]は、望ましい信頼区間になります。

+0

優れた応答をいただきありがとうございました。 Rでは、データ構造全体(そして明示的なモデルさえも)を統計計算機能に渡すことで、同様の方法で実装されているようです。http://www.statmethods.net/advstats/bootstrapping.html – ToddP

+3

ありがとうこれを持ち出す;ちなみに、オーバーフローが発生する可能性が新しくなったように見えるので、左にあるチェックマークを使って良い回答を受け入れると便利だと言えるはずです。 – cge

+0

@ user2269232あなたは答えを受け入れることができます。クリックや2で誰も殺されません.... – rll