階乗関数、二項係数
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
//間瀬茂
SIZE(25){COLOR(red){ 階乗関数、二項係数の計算}}
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
終了行:
//間瀬茂
SIZE(25){COLOR(red){ 階乗関数、二項係数の計算}}
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
ページ名: