Fisher の正確確率検定を Excel を使って計算する

統計学用語辞典

Fisher の正確確率検定を Excel を使って計算する

質問： Excel で，a=63, b=79, c=80, d=64 のようなデータに対して Fisher's exact test の計算をしようとすると，#NUM!と表示されて途中からおかしくなって表示されます。 しかし，../../exact/exact.html の「Fisher's exact test (Extended)」ではちゃんと計算できます。 =FACT( ) を使って計算しており，小さな数字では可能なのですが・・・


問題点は，例数の大きいときの計算がうまく行かないと言うことです。
その原因は，「コンピュータが扱える数の大きさには限界がある」ということです。
記号の定義: 2×2 分割表の 9 個のセルの数値に以下のように記号を割り当てます。

　　　　　　　カテゴリー１　　カテゴリー２　　合計
カテゴリー１　　　　ａ　　　　　　　ｂ　　　　ｅ
カテゴリー２　　　　ｃ　　　　　　　ｄ　　　　ｆ
　　合計　　　　　　ｇ　　　　　　　ｈ　　　　ｎ

教科書や私のページにおいて，生起確率を求めるために，階乗を使う計算式が示されていることがあります。
生起確率 = ( e! × f! × g! × h! ) / ( n! × a! × b! × c! × d! )
しかし，階乗の計算はコンピュータではすぐに（！）オーバーフローしてしまいます。
エクセルでは，170 の階乗 = FACT(170) = 7.2574E+306 が上限です。これだと，合計例数 n が 170 までの分割表なら計算可能と思われるでしょう。
しかし，たとえば，FACT(100) × FACT(100) はオーバーフローします。Fisher の正確確率検定では，FACT の掛算と割算が出てきますので途中の計算結果が大きくなりすぎないように，掛算と割算を交互に行うというのが次の解決法です。
生起確率 = e! / n! × f! / a! × g! / b! × h! / c! / d!
しかし，それでも，限界はすぐに訪れます。
次の解決法は，階乗を使わずに計算する方法です。階乗を使う式は，組み合わせを使う式を展開したものです。組み合わせは，nCi = n! / (n-i)! / i! ですが，COMBIN 関数は途中でオーバーフローが起きないように計算されます。たとえば，n = 171, i = 2 のとき，FACT(171) / FACT(169) / FACT(2) は #NUM! となってエラーになりますが，COMBIN(171,2) = 14535 のように何の問題もなく計算できます。
分割表の生起確率は
生起確率 = COMBIN(e,a) / COMBIN(n,g) × COMBIN(f,c)
で計算できます（COMBIN(e,a) × COMBIN(f,c) / COMBIN(n,g) でもいいですが，途中のオーバーフローの可能性を考えると掛算を先にするのは避けた方がいいです）。
これで，かなり大きな分割表まで取り扱えるようになりますが，それでも限界はあります。
次の工夫は，大きな数値を取り扱うために途中の計算に対数を使うというやりかたです。
まずは，以下のように関数を定義します。

Function logfact(n As Integer) As Double
    Dim i As Long
    logfact = 0
    If n > 0 Then
        For i = 1 To n
            logfact = logfact + Log(i)
        Next i
    End If
End Function

この関数は
n! = 1 × 2 × ... × (n-1) × n
の両辺の自然対数をとって，
log(n!) = log(1)+log(2)+...+log(n-1)+log(n)
の右辺を計算するものです。
n! = EXP(logfact(n)) ですが，まだ EXP は使いません（ここで EXP を使ったら何のメリットもないのです）。
階乗を使う生起確率の式も足し算と引き算になりますので，これらを組み合わすと，
生起確率=EXP(logfact(e)+logfact(f)+logfact(g)+logfact(h)-logfact(n)-logfact(a)-logfact(b)-logfact(c)-logfact(d))
で求めることができます（最後に EXP を使うのがこつ）。
この方法は，log 関数の結果は大きな数字にはならないのと，コンピュータは足し算・引き算は掛算や割算に比べて高速に計算できるという特徴があるので，最後に EXP 関数を使っても元が取れるかもしれません（四則演算の中では，割算は一番コストが高いのです）。
ここでもう少し工夫します。logfact 関数で LOG 関数を何回も使うのはもったいないので，階乗を求める数とその結果の表を前もって計算しておき，途中の計算では logfact 関数を使わずにその表を引けばいいです。表を利用する方法であれば，VBA でプログラムを書かないでも，ワークシートの操作だけで Fisher の正確確率検定を行うこともできます。

定義：
table(0) = 0
for i = 1 to n
	table(i) = table(i-1) + LOG(i)
next i

引用：
LOG(8!) = table(8)

生起確率=EXP(table(e)+table(f)+table(g)+table(h)-table(n)-table(a)-table(b)-table(c)-table(d))

実際には，生起確率はたくさんの分割表について求める必要がありますから，計算量を少しでも少なくするのが望ましいかもしれません。 先の計算式で，table(e)+table(f)+table(g)+table(h)-table(n) の部分はどの分割表においても同じなので，毎回この計算をするのではなくて，最初に constant = table(e)+table(f)+table(g)+table(h)-table(n) としておいて，
生起確率=EXP(constant-table(a)-table(b)-table(c)-table(d))
とすればいいでしょう。
その他の注意としては，変数の取ることのできる値の範囲はその変数の型によって異なるということです。
Excel VBA では，integer は -32,768 〜 32,767 の範囲の整数しかとれません。途中の計算結果も含め，これを越えるような場合には，long を使います。しかし，long も -2,147,483,648 〜 2,147,483,647 の範囲の値しかとれません。また，当然ながらこれらは整数しか扱えないので double を使う必要があります。ちなみにdouble は整数も当然使えます（15 桁くらいの整数も使える）。 は，IEEE 64 ビット (8 バイト) の浮動小数点数の変数です。double は，負の値は -1.79769313486231E308 〜 -4.94065645841247E-324，正の値は 4.94065645841247E-324 〜 1.79769313486232E308 の範囲の値をとります（なお，Single というのもありますが，これを使う必要性はほとんどありません。よほどの理由がない限り使わない方がいいです）。
数値計算は，計算式通りに計算すれば答えが出てくるわけですが，コンピュータを使って計算するときには，計算式通りに計算するだけでは答えが出ないことがあります。
以上を考慮して書いた VBA プログラムがありますので，参考にしてください。