パーミュテーションテストの打ち切り：ベルヌーイ試行尤度積分(4)〜Mathematicaで計算

先日来、こんな記事、こんな記事、こんな記事を書いた。
今、ある多数の仮説検定を行い、マルチプルテスティングの補正をする前の段階で $p_0$ という値が得られているとする。
モンテカルロパーミュテーションにて、最小P値を得るとき、n個の最小P値を得たところ、 $p_0$ 以下のものがk個観測されたとする。
$p_0$ の補正後P値の最尤推定値は $\frac{k}{n}$ であるが、補正後P値の推定値はどのくらいだとみつもるのがよいだろうか。信頼区間はどうだろうか。たとえば、95％信頼区間の上限が0.05を切っていれば、補正後P値が0.05以下であろうことを、95％信頼区間の示す確実度で信じればよいし、95%信頼区間の下限が0.05を越えていれば、補正後P値が0.05以下であることは、95%信頼区間の示す確実度で、ありそうもない、と信じればよい
信頼区間は、『両側』の考え方であるが、ひとまず、『片側』の考え方で、の補正後P値が以下である確率が、95％を越えたら、パーミュテーション試行を打ち切り、逆に、以下である確率が、5%を切ったら、やはりパーミュテーション試行を打ち切るとする。
- たとえば、n=100のとき、 $k>=9$ であれば、 $p_0$ の補正後P値が0.05以下である確率は0.05以下であり、 $k<=1$ であれば、 $p_0$ の補正後P値が0.05以下である確率は0.95以上である。今、補正後P値が0.01であるかどうかを判断基準とすると、 $k>=3$ であれば、補正後P値が0.01以下である確率が0.05以下であるが、たとえ、 $k=0$ であっても、補正後P値が0.01以下である確率は0.95以上にはならない。n=298まで繰り返したときに、 $k=0$ であったときに初めて、補正後P値が0.01以下である確率が0.95以上になる。
これをパーミュテーション試行回数の打ち切り基準とするとすると、その回数基準はこんな数値になる。
今、生起確率pであるような事象について、n=1,2,...,100回、繰り返して実行したところ、そのような事象がk=0,1,2,...,10回、起きたとする。この生起確率が $p \le p_0=\{0.01,0.05\}$ であるような確率が、ガウスの超幾何関数の商で表されるのだが、この超幾何関数の値の算出は、正負の交代のために発散しがち(なことは、こちらと、このサイトの超幾何関数関連の諸記事に詳しい)。そこで、Mathematicaを用いることとする。Mathematicaにも、発散の問題はあるが、手でソースをいじるよりも、広範囲で計算可能である。それを試したのがこちら(第１版)とこちら(第２版)のファイル