![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
階乗関数、二項係数の計算
R 1.9.0 より階乗関数 factorial(n)、その対数値 lfactorial(n) が登場したので、以下の解説はほとんど無意味になりました。 ちなみに factorial(x) = gamma(x+1) がその定義です。
R には階乗関数 n!=1 x 2 x ... x n を計算する関数は無く、しばしば初心者を困惑させるようである。R では階乗はガンマ関数*1を使って次のように計算する。
n! = gamma(n+1) # +1 を忘れないこと！
また、二項係数を計算する関数は choose() である。
choose(n,k)
当然 n, k は整数で 0 <= k <= n である必要があるが、整数でない（それどころか負の）値も平気で受け付けるので注意が必要（これは内部的に gamma() 関数を使っているから）。返り値は常に整数化されるようである。
> choose(3, 1.9) [1] 3
二項係数の自己流定義の例
> nCk <- function(n, k) round( exp( lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1) ) ) > nCk(10, 3) [1] 120
gamma() や choose() はベクトル化されているので、例えば次のような使い方ができる。
> choose(9,0:9) [1] 1 9 36 84 126 126 84 36 9 1 > choose(10:19,9) [1] 10 55 220 715 2002 5005 11440 24310 48620 92378
gamma() や choose() はうっかりすると巨大な数になり、桁溢れを起こし易い。尤度計算等では対数値が求まれば良いので、対数値を与える関数 lgamma()、lchoose() が用意されている。これらは対数値を直接計算するもので、後から対数をとっているわけではない。
lgamma # gamma の自然対数値 lchoose # choose の自然対数値
例えば choose(1000000,1000) がどれくらいの大きさかを知りたければ、次のようにする。つまり、百万人の中から千人を抽出する可能性は 1.507833e+3432 通りある。
> x <- lchoose(1000000,1000)/log(10) # 自然対数値を常用対数値に変換 > x [1] 3432.178 > 10^(x-3432) # 仮数部を計算 [1] 1.507833
