超幾何分布
要因 a, b が独立な場合，以下のような分割表が得られる確率を求める
B not(B) sum
A x o m
not(A) (k-x) o n
sum k o (m+n)
超幾何分布により，求める確率は xCm * xCk-x / m+nCk
R の dhyper の引数はちょっと変で（？），x が生じる確率は dhyper(x, m, n, k)
B not(B) sum
A 2 o 10
not(A) (1) o 5
sum 3 o (15)
> x = 2; m = 10; n = 5; k = 3
> dhyper(x, m, n, k)
[1] 0.494505494505494
> choose(m, x) * choose(n, k-x) / choose(m+n, k)
[1] 0.494505494505495
> exp(lchoose(m, x) + lchoose(n, k-x) - lchoose(m+n, k))
[1] 0.494505494505495
しかし，
> x = 200; m = 2000; n = 5000; k = 600
> dhyper(x, m, n, k)
[1] 0.00104601828356522
> choose(m, x) * choose(n, k-x) / choose(m+n, k)
[1] NaN
> exp(lchoose(m, x) + lchoose(n, k-x) - lchoose(m+n, k))
[1] 0.00104601828356506
となり，普通のやり方では精度が足りない。
Rmpfr パッケージを用いて，単純に計算すると以下のようになる。
library(Rmpfr)
precBits=13
dhyper2 = function(x, m, n, k) {
choose2 = function(a, b) {
prod(mpfr(1:a, precBits)) / prod(mpfr(1:b, precBits)) / prod(mpfr(1:(a-b), precBits))
}
as.numeric(choose2(m, x) * choose2(n, k-x) / choose2(m+n, k))
}
> x = 200; m = 2000; n = 5000; k = 600
> dhyper(x, m, n, k)
[1] 0.00104601828356522
> dhyper2(x, m, n, k)
[1] 0.00104601828356522
dhyper2 の実行時間は 0.15 秒
