fc2ブログ
回帰係数の一様性の検定
マニア度があがりすぎて、誰も興味ないかもしれませんが、ちょっと仕事で簡単なスクリプトを作ったので、貼っておきます。
共分散分析をするときに、入力データの各群の回帰係数が等しいかどうかを先に検定します。これで等しくないとその先には進めないのですが、逆にこの部分だけを使っていろいろ物事を言うこともできるかなというのもあり、実装してみました。
武藤眞介著「統計解析ハンドブック」朝倉書店p.p.272-275を参考にしています。273ページの式(4)に間違いがあるので注意してください。(初版6刷)
スクリプトでは修正済みです。いやー、scipy大変便利です。

from scipy.stats.distributions import f

# input data is list of below
# data1 [[x,,,],[y,,,]]
# data2 [[x,,,],[y,,,]]
#....
# return (f_v , p-value )

def test_regcf_uni( data ):
 n = 0.0
 for (i,v) in enumerate(data):
  if len(v[0]) <> len(v[1]):
   raise ArithmeticError, str(i) + 'th data length inconsistency'
  n += float(len(v[0]))
 k = float(len(data))
 u = []
 v = []
 w = []
 r = []
 for (i,_d) in enumerate(data):
  u.append( sum( float(_x)*_x for _x in _d[0]) - (sum(_d[0])*sum(_d[0]) / float(len(_d[0]))) )
  v.append( sum( float(_y)*_y for _y in _d[1]) - (sum(_d[1])*sum(_d[1]) / float(len(_d[1]))) )
  w.append( sum( float(_x)*_y for _x,_y in zip(_d[0],_d[1])) - (sum(_d[0])*sum(_d[1]) / float(len(_d[0]))) )
  r.append( v[i] - w[i]*w[i] / u[i] )
 SR = sum(r)
 SE = sum(v) - sum(w)*sum(w) / sum(u)
 f_v = ( (SE-SR)/(k-1.0) )/( SR/(n-2.0*k) )

 return (f_v, f.sf( f_v , k-1 , n-2*k ).tolist())

sample_data = [ [[4,2,1,3,5],[180,150,135,175,200]],\
       [[3,5,6,4,7],[155,190,205,165,210]],\
       [[3,6,2,5,4],[190,240,175,215,205]] ]

print test_regcf_uni( sample_data )

改行されてしまっているところがありますが、¥が入っていなければ、実際は1行です。
そのまま、切り貼りでお使いいただいてかまいませんが、もしよく分からなかったら、ご一報ください。
スポンサーサイト



テーマ:数学 - ジャンル:学問・文化・芸術

【2008/07/01 16:24】 | 数学 | トラックバック(0) | コメント(0) | page top↑
Pythonでの繰り返し操作
別にネットを探せばどこにでもある情報なのであまり価値はないのですが、こうして至る所にPythonがイイって書いてあれば、Pythonが日本でもっと普及するかなと思い、いくつか便利な点をまとめておきます。

ある変数の値の範囲を調べる時、Pythonでは、
In [4]: g = 2
In [5]: 1 < g < 5
Out[5]: True
のようにできます。書くときのストレスがないです。

In [1]: l = [1,2,3,4]
In [2]: r = [2,2,2,2]
とすると、
In [9]: [ x*y for x in l for y in r]
Out[9]: [2, 2, 2, 2, 4, 4, 4, 4, 6, 6, 6, 6, 8, 8, 8, 8]
これは、lとrのループを総当たりでまわすことに相当します。

In [10]: [ x*y for x,y in zip(l,r)]
Out[10]: [2, 4, 6, 8]

こう書くと、2つのリストを対応づけてくれます。

In [12]: for (i,x) in enumerate(r):
....: print i,x
0 2
1 2
2 2
3 2

とすると、リストの要素に番号をふってくれます。iteraterを利用した表記は便利だけど、操作のためにインデックスが欲しいときに便利です。引数はsetでも大丈夫。

ちなみに、リスト内包表記は、英語ではlist comprehension
こういう表記の揺らぎを考慮した検索エンジンそろそろ出ても良さそうな気が。

テーマ:プログラミング - ジャンル:コンピュータ

【2008/07/01 13:34】 | Python | トラックバック(0) | コメント(0) | page top↑
MacOS10.5.3にscipy0.6.0をインストールする
numpy-1.1.0を入れたときは、すんなり入ったような気がしていたので、簡単に入るかなと思ったのですが、ソースをコンパイルしてインストールしたら、なぜかうまく動かない。
という訳で、ちょっと探すと、MacOSのみなさまのためのページが。
http://www.scipy.org/Installing_SciPy/Mac_OS_X
お、これで言われたとおりにすればいいとおもったんですが、なぜかfftw-3.1.2のmakeに失敗します。configureのあと、makeで即エラー。悲しい感じですが、
sudo port install fftw-3
で入りました。
ちなみに、Mac Portsは好きなのですが、たとえばscipyをインストールしようとすると、なぜかシステムのPythonを見つけられずに、Pythonのコンパイルから始めたりして、最初から使わない奴は無視か?と思わせるところが残念なのですが、勘違いかな。ちなみに、python.orgの2.5.2とAppleの2.5.1の両方入っています。
ひとつ気をつけないといけないのは、scipyのインストールが完了したあと、pythonを起動して、import scipyするのは、scipyのソースディレクトリを抜けてから。前述のページに書いてありますが、やりがちです。

テーマ:フリーソフト - ジャンル:コンピュータ

【2008/07/01 00:54】 | Python | トラックバック(0) | コメント(0) | page top↑
| ホーム |