マニア度があがりすぎて、誰も興味ないかもしれませんが、ちょっと仕事で簡単なスクリプトを作ったので、貼っておきます。
共分散分析をするときに、入力データの各群の回帰係数が等しいかどうかを先に検定します。これで等しくないとその先には進めないのですが、逆にこの部分だけを使っていろいろ物事を言うこともできるかなというのもあり、実装してみました。 武藤眞介著「統計解析ハンドブック」朝倉書店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行です。 そのまま、切り貼りでお使いいただいてかまいませんが、もしよく分からなかったら、ご一報ください。 スポンサーサイト
|
|
| ホーム |
|