階乗関数、二項係数の計算

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

*1 注: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)=√π。

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1332d)