一般化正規分布をStanで定義してMCMCしてみる

誰得なのかよくわからない記事を書きます。※最後に追記あります

独立成分分析を調べていたら,ラプラス分布と正規分布の両方を含んだ「一般化正規分布」というのが出てきたので,調べてみました。

一般化正規分布は,平均μと,尺度パラメータα,形状パラメータβの3つのパラメータがあります。この分布の特徴は,左右対称の分布で,尖度を形状パラメータβで変えることができるというのが特徴です。βが1のときはラプラス分布,βが2のときは正規分布,βが∞のときはμ-αとμ+αの範囲の一様分布となります。

ついでに,ラプラス分布というのは,平均値あたりがすごいとんがっている分布です。

ggd2

一般化正規分布。楽しそうですね。

しかし,Stanにこの分布は入っていません。なので前の記事の方法を使って,StanでMCMCできるようにしてみました。

まず,一般化正規分布の確率密度関数は次のようになります。

ggd1

Γ()はガンマ関数です。

まず,Stanでfunctions{}ブロックに一般化正規分布の上の確率密度関数を書きます。関数名はggd_log()です。名前の通り,戻り値はlogしたものにします。また,Stanでは絶対値はabs()は推奨されて無くて,fabs()のほうがよいようです。イマイチ違いはわかってません。

あとは,この分布を使ってパラメータを推定するようなStanコードを書けばOKです。もちろん,x ~ ggd()というようにサンプリング関数として使えます。

下のコードでは,一般化正規分布の標準偏差をgenerated quantities{}で吐き出すコードも加えています。

それでは,正規乱数から生成したデータで推定してみましょう。なお,上のコードはfunction2.stanという名前で保存しています。

平均5,標準偏差2の正規乱数を1000個作り,推定してみました。

結果は次のようになりました。

ほぼmu=5.03,σ=1.99とほぼ母数通りとなりました。またβは1.96となり,2に近い値が推定されました。

次にラプラス分布を生成してみます。

ラプラス分布をして推定してみます。パラメータは平均μと尺度パラメータφです。φは一般化正規分布のαと同じです(たぶん)。

推定した結果が以下です。

μもαも母数に近い値となっています。βは1.02と,ラプラス分布の値とほぼ一致しています。

最後に一様分布です。0から2の範囲の一様分布から1000個生成しました。

結果は,以下です。

μが1でαが1なので,0~2の一様分布として推定できているのがわかります。また一様分布の場合βは∞となりますが,今回は372と推定されています。十分大きい値なのでだいたいOKという感じでしょうか。

このように,一般化正規分布は,ラプラス~正規~一様分布という3つの分布がβという形状パラメータで連続的に表現される楽しい分布でしたとさ。

 追記

beroberoさんから,確率密度は最後のlogとるんじゃなくて,最初からできるだけlog使って計算したほうが速いよ!というご指摘をもらいました。

 というわけでご指摘にようにコードを修正してみました。

すると,改善前のコードだと,13秒かかってたのが8.5秒と,かなりスピードアップしました。というわけで,こっちのほうがオススメです。

This entry was posted in 心理統計学, Stan. Bookmark the permalink.