//間瀬茂 SIZE(25){COLOR(red){ 階乗関数、二項係数の計算}} COLOR(red){R 1.9.0 より階乗関数 factorila(n)、その対数値 lfactorial(n) が登場したので、以下の解説はほとんど無意味になりました。} ちなみに factorial(x) = gamma(x+1) がその定義です。 COLOR(red){R 1.9.0 より階乗関数 factorial(n)、その対数値 lfactorial(n) が登場したので、以下の解説はほとんど無意味になりました。} ちなみに factorial(x) = gamma(x+1) がその定義です。 ~ ~ R には階乗関数 n!=1 x 2 x ... x n を計算する関数は無く、しばしば初心者を困惑させるようである。R では階乗はガンマ関数((注:x>0 ならば Γ(x) は t^(x-1)exp(-t) を 0 から ∞ まで積分したものと定義される。解析接続により Γ(x) は x=0,-1,-2,... を除くすべての実数に対して拡張定義できる。R のΓ関数はこの拡張された Γ 関数である。基本的な性質は 1) x, x+1 が共に正でない整数でなければ Γ(x+1)=xΓ(x), 2) n が自然数ならΓ(n+1)=n!, 3) Γ(1)=1, 4) Γ(1/2)=√π。))を使って次のように計算する。 n! = gamma(n+1) # +1 を忘れないこと! また、二項係数を計算する関数は SIZE(20){choose()} である。 choose(n,k) 当然 n, k は整数で 0 <= k <= n である必要があるが、整数でない(それどころか負の)値も平気で受け付けるので注意が必要(これは内部的に SIZE(20){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 SIZE(20){gamma()} や SIZE(20){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 SIZE(20){gamma()} や SIZE(20){choose()} はうっかりすると巨大な数になり、桁溢れを起こし易い。尤度計算等では対数値が求まれば良いので、対数値を与える関数 SIZE(20){lgamma()}、SIZE(20){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