//間瀬茂
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


トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS