No.00505　P値　　【まな】　2006/06/29(Thu) 10:40

はじめまして。現在自分でT検定，ANOVAのプログラムを作成しようと考えています。T，自由度までは計算式をみつけ出せるようになりましたが，そこからP値を計算するにあたらり，
分布表をつかわずに計算式で出したいのですが，このような計算式を教えてはいただけないでしょうか？何卒お願いします。
No.00506　Re: P値　　【青木繁伸】　06/06/29(Thu) 11:33

どのような言語でプログラムしているのでしょうか？
場合によっては，ライブラリが用意されていると思いますので。
No.00507　Re: P値　　【まな】　06/06/29(Thu) 12:17

ご連絡ありがとうございます。dellphiですが，計算式自体も教えていただければ幸いです。
遺伝子の 解析を行っており，サンプル数や内容により，student t-test，Welchit-test，oneway-ANOVA,two-way-ANOVAを行いたいと考えています。さらに，数万のデータを扱う ので多重検定も後々考えています。お手数ですが，ご教授のほどよろしくお願いいたします。
No.00508　Re: P値　　【ひの】　06/06/29(Thu) 13:23

Pascalでの数値計算ライブラリはネットで探すといろいろあります。私が使ったの（最近はプログラミングしていないので過去形ですが）は，TP Math というライブラリです。Turbo Pascalの時代からあり，現在もアクティブなページなのでバグは枯れているのではないかと思います。マルカート法のプログラミングに挫折して(^^;)このライブラリのお世話になりました。主な分布関数の計算も含まれていたと思います。URLは以下。

http://www.unilim.fr/pages_perso/jean.debord/tpmath/tpmath.htm

その他リンク集としては，

http://www-rab.larc.nasa.gov/nmp/nmpIndex.htm
No.00509　Re: P値　　【青木繁伸】　06/06/29(Thu) 15:15

アルゴリズムということでしたら，いろんな言語で書かれていても，翻訳できればいいわけで
http://aoki2.si.gunma-u.ac.jp/Hanasi/Algo/getprob.html
などもありますが，個々のアルゴリズムがどこに根拠があったのかはもう不明になりました。

Numerical recipes in C
などは，権威があるので，そのあたりを見るのもよいかと
No.00559　Re: P値　　【青木繁伸】　06/07/05(Wed) 22:58

強制的にスレッドを維持します

-------------------------------------------------------------------------------

P値　投稿者：まな 投稿日：2006/07/05(Wed) 14:59 No.554

先日ご教授いただいたとおりプログラムから抜粋した計算式でP値の計算をしました。
Student-t-testの式より求めたT=2.365，自由度＝8をP値計算式にはめ込みました。その結果，私の計算値ではP＝0.045456355となりました。しかし，エクセルから算出したPの値は0.045600927，その他別のフリーソフトで行ったものは0.04563694893441，または0.04560092748824となります。ほかは0.00456まで同じですが私の場合は0.0454となります。どこか間違っているのでしょうか？また，完全に同じP値はなかったのですが，どれが正しいとか，どの部分が違うため値が変わってくるなど教えていただけないでしょうか？
No.00560　Re: P値　　【青木繁伸】　06/07/05(Wed) 23:19

スレッドを維持するように，先の発言への返信という形で投稿して欲しかったですね。

貴方が使った近似式がどれであるかと，貴方が実際に書いたプログラムが提示されていないので，どこが間違えているか，指摘のしようがありませんね。

近似式は，その精度もいろいろありますので，全く同じ答えにならない場合がありますが，，，ちょっと違いすぎでしょう

R だと
> options(digits=20)
> pt(2.365, 8, lower=FALSE)*2
[1] 0.04560092748838392
ということになるでしょう。

もう一つの可能性は「Student-t-testの式より求めた t=2.365」って，正確にその値になったのでしょうか？
貴方が表示したときに 2.365 って表示されただけで，内部的には異なる値なのでは？
貴方が書いたプログラムに t=2.365，df=8 を与えたらどうなるんですか？
No.00561　Re: P値　　【青木繁伸】　06/07/05(Wed) 23:20

私が紹介した
http://aoki2.si.gunma-u.ac.jp/Hanasi/Algo/getprob.html
は，以下のようになりますね。
aoki@aoki1 [9] > cat t.awk
function abs(x)
{
   return x < 0 ? -x : x
}
function txp(t, df,    i, t1, t2, p, w, m2pi)
{
   m2pi = 0.636619772367581343076
   t1 = abs(t)/sqrt(df)
   t2 = 1.0/(1.0+t1*t1)
   if ((df%2) == 0) {
       w = t1*sqrt(t2)
       p = 1.0-w
       for (i = 2; i <= df-2; i += 2) {
           p -= w *= t2*(i-1)/i
       }
   }
   else {
       p = 1.0-m2pi*atan2(t1, 1)
       if (df >= 3) {
           w = m2pi*t1*t2
           p -= w
           for (i = 3; i <= df-2; i += 2) {
               p -= w *= t2*(i-1)/i
           }
       }
   }
   return (p < 0.0 && abs(p) < 1e-10) ? 0.0 : p
}
BEGIN {
        printf("%.16g\n", txp(2.365, 8))
}

aoki@aoki1 [10] > gawk -f t.awk
0.04560092748838412
No.00566　Re: P値　　【まな】　06/07/06(Thu) 10:41

青木先生，お返事ありがとうございます。
私が行った計算式ですが，以下のようなものです。


Student’s t-testの式より
　　　　　　　T=2.365
df=8
以下，p値計算の式に当てはめると，
　　　　　　　t1=2.365/√8=0.836
t2 = 1.0/(1.0+(0.836)2)=0.589
w1=0.836*√(0.589)=0.641
p1=1.0−0.641=0.359
次にi=2となるので，
　　　　　　　w2=0.641*0.589*((2−1)/2)=0.189
p2=0.359−0.189=0.170
この後，iが2づつ増えていくと，df−2=6になるまでにあと2回上の計算を繰り返すことになる。
　　　　　　　i=4 w3=0.189*0.589*(3/4)=0.0835
p3=0.170−0.0835=0.0863
i=6 w4=0.0835*0.589*(5/6)=0.0410
p4=0.0863−0.0410=0.0453
p4の値(0.0453)が最終的なp値となる。(Excelでは0.0456)


以上が私の行った計算です。

T=2.365, df=8
の2つの値をExcelのTDIST関数や他の解析ソフトに代入するとp=0.0456程度になりました。
また，上記の計算式の有効数字を増やしても0.0454ぐらいにしかなりません。

上記の計算式のどこかに間違いがあるのでしょうか？
ご教授頂ければ幸いです。
No.00567　Re: P値　　【青木繁伸】　06/07/06(Thu) 10:59

> p値計算の式に当てはめると，
> t1=2.365/√8=0.836
> t2 = 1.0/(1.0+(0.836)2)=0.589
> w1=0.836*√(0.589)=0.641
> p1=1.0−0.641=0.359

え？このように手計算（コンピュータは使うものの）したのですか？
t1 が 0.836 で良いわけがありません。途中の計算で小数点以下3桁しか使わなければ結果が小数点以下3桁より正確な答えは出ないでしょう？
コンピュータは内部に有効桁数は16桁ほどはあるんですよ。
t1 ≒ 0.8361537687530924
手計算して確かめるなら，有効桁数を十分に大きく取っておかないと。。。

R で行った計算例を示しましょう
R では，結果が逐次的に表示できるので，こういうときには便利ですが。
（ ）でくくると，その中に示された計算式（代入）の結果が表示されるんですね。
表示されている数値は，内部で保持されている数値です。

> options(digits=22)
> (t=2.365)
[1] 2.3650000000000002
> (df=8)
[1] 8
> (t1 = abs(t)/sqrt(df))
[1] 0.8361537687530924
> (t2 = 1.0/(1.0+t1*t1))
[1] 0.588528476502081
> (w = t1*sqrt(t2))
[1] 0.641460461367588
> (p = 1.0-w)
[1] 0.358539538632412
> (i = 2)
[1] 2
> (w = w*t2*(i-1)/i)
[1] 0.18875887403249428
> (p = p-w)
[1] 0.16978066459991767
> (i = 4)
[1] 4
> (w = w*t2*(i-1)/i)
[1] 0.08331747942044405
> (p = p-w)
[1] 0.08646318517947361
> (i = 6)
[1] 6
> (w = w*t2*(i-1)/i)
[1] 0.040862257691089518
> (p = p-w)
[1] 0.0456009274883841

毎回の演算結果を3桁でやれば結果に大きな誤差が含まれます。
コンピュータを使うということは，途中結果をどの程度の桁数まで使うかを考える必要なんかは本来はないのですから。
No.00582　Re: P値　　【まな】　06/07/06(Thu) 21:30

お返事ありがとうございます。

どうやら私の計算方法に問題があったようです。
コンピュータを使ってはいましたが，有効数字については深く考えずに計算を行っていました。
お忙しい中，どうもお手数をおかけして申し訳ありませんでした。


再び質問で恐縮ですが，小数自由度を持つ Welchのt-test のp値については

http://aoki2.si.gunma-u.ac.jp/Hanasi/Algo/getprob.html

のt分布における計算式は適用できないのでしょうか？
Welchの式から求めたtと自由度を代入して計算を行ったところ，ExcelのTTEST関数から求めたp値とは前以上の違いがありました。
今回は有効数字も考慮に入れ，コンピュータを使用したので問題はないはずです。
小数自由度の場合は，他の手法が必要なのでしょうか？
ご教授頂ければ幸いです。
No.00583　Re: P値　　【青木繁伸】　06/07/06(Thu) 22:12

小数自由度を認めるかどうかも，近似公式によるわけですよ。
件の近似公式は，小数自由度は考慮していません。（たぶん）

小数自由度も許容する近似公式については，私が補足しているものはありませんが，Excel で計算することについては，投稿頂いた記事があり
http://aoki2.si.gunma-u.ac.jp/Hanasi/StatTalk/arc-s001.html
を参照して頂けると分かると思います。
これを読めば分かるように，ベータ関数とガンマ関数を用意してしかるべく計算すれば，小数自由度を持つカイ二乗分布，t分布，F分布であろうと，ちゃんと計算できます。

なお，この件からもよく分かるように，既にあるプラットフォームで確立している計算方式（大げさに言えば，技術）を，苦労して別のプラットフォームに移植するのは意義のあることでもありますが，既に十分に確立・検証されており，しかも付加的能力が十二分にあるプラットフォームがあるなら，あえてそれ以外のプラットフォームを構築する意義がどれほどあるか。

簡単に構築できるなら良いですが，苦労に苦労を重ねて，なおかつ元の何分の一しか実装できないならむなしいではないですか。　。。。あ，ちょっと，余計なことをいいました。
No.00585　Re: P値　　【まな】　06/07/07(Fri) 00:05

お返事ありがとうございます。

どうやら自力で計算することにこだわり過ぎていたようです。
今後はご助言頂いた通り，すでに確立された方法については，それらを最大限活用するようにしたいと思います。

お付き合い頂いて誠にありがとうございました。

● 「統計学関連なんでもあり」の過去ログ--- 038 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る