ベイズマックスへようこそ。
このブログでは,心理学と統計が大好きなmutopsyが,
統計に関するあれやこれやを自由気ままに投稿します。

2017年12月30日

Stanでポン・デ・リングを作ってみた

 こんにちは。mutopsyです。この記事は,Stan Advent Calendar 2017の25日目のエントリー記事です(実際にこの記事を書いたのは12月30日ですが)。

 みなさんはドーナツは何が好きですか?僕はフレンチクルーラーとエンゼルクリームがお気に入りですが,たまーにポン・デ・リングなんかも食べたくなります。このブログを訪れてくださった皆様の中には「Stanでポン・デ・リングが作れたら嬉しいのになあ……」とお考えになられた経験のある方もおられるでしょう。本日,なんとなくTogetterまとめを見ていたら,こんなまとめを見つけました。

「おっぱい関数選手権」の名古屋代表の名大生さん、1つの数式だけでポン・デ・リングを表示することに成功「天才の所業」

 このまとめは,ポン・デ・リングの形状(三次元と二次元の両方)がたった1つの数式で表現可能であることを発見した@CHARTMANqさんのツイートをまとめたものです。このまとめを読んで甚く感動した僕は,さっそくStanで実装してみることにしました。Stanを乱数生成器として用いて,二次元のポン・デ・リング分布からのMCMCサンプルを発生させるという試みです。三次元でもできますが,とりあえず二次元ポン・デ・リング不等式を用いてみました。

 二次元ポン・デ・リング不等式は,極座標を用いて(r-2|cos4θ|-10)(r+|cos4θ|-6)≤0と表せます(上記まとめの@CHARTMANqさんのツイートをご参照ください)。極座標というのは,原点からの距離 (r) と角度 (θ) を使って空間上(ここでは平面上)の点の位置を表す座標系です。この不等式を満たすrとθの領域が 二次元のポン・デ・リングを表します。

 二次元ポン・デ・リング分布からの乱数(MCMCサンプル)を発生させるStanコードの例は以下の通りです。

Stanコード01: 二次元のポン・デ・リング分布から乱数を発生させるコード
parameters{
  real<lower=0> r;  //極座標のr
  real<lower=0,upper=2*pi()> theta;  //極座標のθ(0°〜360°の範囲)
}

transformed parameters{
  real<upper=0> D;  //上限を指定することで不等式を表現
  D = (r - 2*fabs(cos(4*theta)) - 10)*(r + fabs(cos(4*theta))-6); //二次元ポン・デ・リング(不)等式
}

generated quantities{
  real x;
  real y;
  x = cos(theta) * r;  //直交座標のxを計算
  y = sin(theta) * r;  //直交座標のyを計算
}

 今回はStanを乱数生成のためだけに用いているのでdataブロックは必要ありません。また,事前分布として一様分布を用いるため,modelブロックも不要です。Stanに慣れていると変な感じがするかもしれません。rとθの事前分布(一様分布)のうち,ポン・デ・リング不等式を満たすパラメータのみをMCMCサンプルとして生成することにより,二次元ポン・デ・リング一様分布からの乱数を発生させます。このモデルを実行するためのRコードは以下の通りです。

Rコード01: 上記モデルを走らせるためのRコード
library(rstan)
data <- list(Dummy = 1) #ダミーデータを作る(実際には使わない)
stanmodel <- stan_model(file="stancode/Ponde_2d.stan") #上記Stanコードを読み込む
fit <- sampling(
  stanmodel,data=data,seed=123,
  chains=4,
  iter=10000
)

 上記のコードではダミーデータを渡していますが,実際にはこのデータに意味はありません。Stanコードを走らせてみたところ,Rhatは全て1.01未満となり問題なく収束しました。各パラメータ (r, θ, x, y) の事後分布は以下のようになります。

ponde_marginal.png

 周辺分布だけを見ても,なんとなくポン・デ・リングっぽくありませんか。θの分布の凹凸やxとyの二峰分布具合はまさにポン・デ・リングですよね。ではいよいよxとyの同時分布を見てみましょう。

pon.png

 どう見てもポン・デ・リングです。本当にありがとうございました。ぜひ皆さんもStanでポン・デ・リングを作ってみてください。もちろん数式が分かればポン・デ・リング以外のものも作れます。お試しあれ。それでは皆様,よいお年を!

 (「別にStanじゃなくてもよくね?」とつっこんではいけません。)

posted by mutopsy at 21:15 | 雑記

2017年12月07日

統計ブログを作りました。

こんにちは。mutopsyこと武藤拓之です。
大阪大学で認知心理学の研究をしています。
心理学と統計が好きな,どこにでもいる普通のベイジアンです。

勢いあまって統計ブログというものを作ってしまいました。
RやStanを使って統計に関する記事を書いたり,あるいは書かなかったりする予定です。

このブログの名前はベイズマックス。あなたの知的好奇心をくすぐります。さういふブログに私はしたい。
何卒,よろしくお願いいたします。
posted by mutopsy at 16:08 | 雑記