概要
正規分布の場合、x^nを積分しなくても、母関数をn階微分しなくても、原点周りの高次モーメントが順次求まります。
3次モーメント・4次モーメントのみ前回解説済み。
本題
前回、正規分布の3次モーメントE[x^3]は
のように求められることを紹介しました。これは漸化式だったりします。
Rのpracma/quadprogというライブラリを使うと、多項式微分(POLYnomial DERivative)ができ、出力を係数ベクトルとして得られますので、この漸化式的なものを解くことができます。
xはμのn次に対する係数ベクトルなので、c(x, 0)というのが、μを掛け算する操作に相当します。
1 2 3 4 5 6 7 8 9 10 11 |
install.packages(pracma) install.packages(quadprog) library(pracma) f <- function(x) c(0, 0, polyder(x)) + c(x, 0) x <- c(1, 0) # 初期値 E[x] = μ(+0) for (i in 1:10) { print(x) x<-f(x) } |
これだけで完了です。実行結果を記します。
1 2 3 4 5 6 7 8 9 10 |
[1] 1 0 [1] 1 0 1 [1] 1 0 3 0 [1] 1 0 6 0 3 [1] 1 0 10 0 15 0 [1] 1 0 15 0 45 0 15 [1] 1 0 21 0 105 0 105 0 [1] 1 0 28 0 210 0 420 0 105 [1] 1 0 36 0 378 0 1260 0 945 0 [1] 1 0 45 0 630 0 3150 0 4725 0 945 |
この出力は以下を意味しています。
・・・もちろん、モーメント母関数を10回微分していただいても結構でございます。
ピンバック: 正規分布の3次モーメントの楽な計算法(確率論の裏ワザ) | teqニカルブログ