正規分布のn次(高次)モーメントを生成するプログラム(続・確率論の裏ワザ)

概要

正規分布の場合、x^nを積分しなくても、母関数をn階微分しなくても、原点周りの高次モーメントが順次求まります。
3次モーメント・4次モーメントのみ前回解説済み。

本題

前回、正規分布の3次モーメントE[x^3]は

E[x^3] = \sigma^2 \partial_\mu E[x^2] + \mu (\mu^2 + \sigma^2) \\= \sigma^2 (2 \mu) + \mu^3 + \mu \sigma^2 \\= \mu^3 + 3 \mu \sigma^2

のように求められることを紹介しました。これは漸化式だったりします。

E[x^{n+1}] = \sigma^2 \partial_\mu E[x^n] + \mu E[x^n]

Rのpracma/quadprogというライブラリを使うと、多項式微分(POLYnomial DERivative)ができ、出力を係数ベクトルとして得られますので、この漸化式的なものを解くことができます。

xはμのn次に対する係数ベクトルなので、c(x, 0)というのが、μを掛け算する操作に相当します。

これだけで完了です。実行結果を記します。

この出力は以下を意味しています。

E[x] = \mu \\E[x^2] = \mu^2 + \sigma^2 \\E[x^3] = \mu^3 + 3 \mu \sigma^2 \\E[x^4] = \mu^4 + 6 \mu^2 \sigma^2 + 3 \sigma^4 \\E[x^5] = \mu^5 + 10 \mu^3 \sigma^2 + 15 \mu \sigma^4 \\E[x^6] = \mu^6 + 15 \mu^4 \sigma^2 + 45 \mu^2 \sigma^4 + 15 \sigma^6 \\E[x^7] = \mu^7 + 21 \mu^5 \sigma^2 + 105 \mu^3 \sigma^4 + 105 \mu \sigma^6 \\E[x^8] = \mu^8 + 28 \mu^6 \sigma^2 + 210 \mu^4 \sigma^4 + 420 \mu^2 \sigma^6 + 105 \sigma^8 \\E[x^9] = \mu^9 + 36 \mu^7 \sigma^2 + 378 \mu^5 \sigma^4 + 1260 \mu^3 \sigma^6 + 945 \mu \sigma^8 \\E[x^{10}] = \mu^{10} + 45 \mu^8 \sigma^2 + 630 \mu^6 \sigma^4 + 3150 \mu^4 \sigma^6 + 4725 \mu^2 \sigma^8 + 945 \sigma^{10}

・・・もちろん、モーメント母関数を10回微分していただいても結構でございます。

「正規分布のn次(高次)モーメントを生成するプログラム(続・確率論の裏ワザ)」への1件のフィードバック

  1. ピンバック: 正規分布の3次モーメントの楽な計算法(確率論の裏ワザ) | teqニカルブログ

コメントは受け付けていません。