Kosugitti's BLOG

アンドロイドは正規分布の夢を見るか

この画面は、簡易表示です

stan

習作;混合正規分布モデル(Gaussian Mixture Model) by rstan

Stanのマニュアルに載っている簡単なモデルなんだけど,最後の「各要素の所属確率」を算出するのになぜか手間取った。softmax関数を使うときの型の問題でした。ちぇ。

Stanコードはこの通り。generated quantitiesのところでsoftmax関数を使います。

data{
  int<lower =2> K; #number of clusters
  int<lower =1> N; #number of observations
  real X[N]; #observed data
}

parameters{
  vector[K] mu;
  vector<lower =0,upper=10>[K] sig2;
  simplex[K] theta;
}

transformed parameters{
  vector[K] ps[N];
  for(n in 1:N){
    for(k in 1:K){
      ps[n,k] = log(theta[k])+normal_lpdf(X[n]|mu[k],sig2[k]);
    }
  }
  
}

model{
  sig2 ~ cauchy(0,2.5);
  mu ~ normal(0,100);
  for(n in 1:N){
    target += log_sum_exp(ps[n]);
  }
}

generated quantities{
  simplex[K] u[N]; #class membership probability
  for(n in 1:N){
    u[n] = softmax(ps[n]);
  }
}

キックするコードは次の通り。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# 1 dimensional Mixed Gauss
#  K個の正規分布を混ぜる

K < - 3
mean1 <- 5
mean2 <- -5
mean3 <- 0
sd1 <- 1
sd2 <- 2
sd3 <- 3

N1 <- 500
N2 <- 300
N3 <- 200
N <- N1+N2+N3
X <-c(rnorm(N1,mean1,sd1),rnorm(N2,mean2,sd2),rnorm(N3,mean3,sd3))

### MCMC!
stanmodel <- stan_model("develop/gmm_single.stan",model_name="GMM")
standata <- list(K=K,N=N,X=X)
max.iter <- 3000
burn.in <- 1000
itr <- max.iter - burn.in
C <- 4


# vbは速いしラベルスイッチングとかないから嬉しい
fit_vb <- vb(stanmodel,data=standata)
print(fit_vb,pars=c("mu","sig2","theta"))



# 一つならラベルスイッチングは起きない
fit_sp <- sampling(stanmodel,data=standata,chain=1,iter=max.iter) 
# Rhatも優秀
print(fit_sp,pars=c("mu","sig2","theta"))
print(fit_sp,pars=c("u"))


# チェインを複数にするとラベルスイッチングが起きる 
fit_sp <- sampling(stanmodel,data=standata,chain=C,iter=max.iter,warmup=burn.in)
# Rhatが悪くなってもくじけちゃいけない
print(fit_sp,pars=c("mu","sig2","theta"))
# 目で見ると綺麗なラベルスイッチングが確認できるよ。
traceplot(fit_sp,pars="mu")


塾合宿のおもひで

お盆休みを満喫しております。

 

お休みに入る前,8/8-8/11の二泊三日で魁!ベイズ塾夏合宿2016に参加してきました。今年から新塾生も加入し,総勢19名というソコソコの集団サイズ。残念ながら,最近の若手は魁!男塾を知らないようで,そこかしこにある「塾のオマージュ」が伝わっていないようです。

二日間にわたって,StanとRを使いながら,ベイズ推定に色々な角度からアプローチ。日中は講義スタイル,夕食後はLTで自主的に研究発表をするというスタイルだったのだけど,さすが男気あふれる男塾,みんな自分で面白いことをやってますなあ。私はゼミ生のモデルの紹介と,状態空間モデルの話を披露しました。個人的にはアニメの公式ツイッターアカウントの,フォロアー数を自動的にとってこれるAPIの存在を知ったことがでかかったかな。いや,もちろん一番は国里先生の機械学習の話でしたけどね。

今後の活動について,どうしましょうかという話が最後にありましたけど,基本的には塾生のスキルアップ,塾からの情報発信だと思うので,塾長に何か面白いことを思いついてもらうかなぁ。最近はいろんなところでベイズの理解者が増えてきたとは思うので,目新しいことをするのが難しくなってきたかなとは思うね。



毒くらわば皿まで

rstan 2.10が来てますね。

2.9とちがって,式の代入が< -から=にすることが推奨されているようで,そうなっちゃうと今までのコードも全部書き換えないといけない。
とはいえこれから先はこうなるってんなら,早めに対応しておこうか,と母艦などのパッケージを思い切ってアップデート!

本艦のMacはすんなりアップデートでき,サンプルとして昔ながらの8schoolsを試したら,無事に走りました。もっともエラーとして,


DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal): assignment operator < - deprecated in the Stan language; use = instead.

って言われちゃいましたけどね。あー,やっぱりかー(当たり前だ)と。

で,書き直しなんだけど,よく見たら


theta[j] = mu + tau * eta[j];

だけでなく,モデルセクションも


model {
target += normal_lpdf(eta | 0, 1);
target += normal_lpdf(y | theta, sigma);
}

って変わってるのね。もとは


model {
eta ~ normal(0, 1);
y ~ normal(theta, sigma);
}

だったんですよ。

ともかくこれで走りました。コンパイルの時にちょっと警告が出たけどね。
これからは,コンパイル時のこの警告を潰していかないといけません。うう,やっと関数が出来上がったと思ったら・・・(苦笑)

[crayon-5bc7040c5f06f648500891/]

 

ついでに,とMacの上の仮想Windowsは,っとみたら,Rが3.2で止まってた。R3.3.1にあげて,rstanもアップデートする。

ここでRtoolsのバージョン違いがあってうまくいかなかった。最新版Rtools34はダメみたいで,Rtools33にしか対応してないみたい。でそこからが大変。コンパイル系のエラーが出まくって,rstanがうまく入りません!

仕方がないので,RもRtools33もRstudioもぜんぶアンインストール&再インストールしてから,install.packages("rstan")した。コンパイル時になんか警告が出たっぽいけど,これでやっとインストールできました。
ちなみにrstudioapiとかcolorなんとかってパッケージが入らないみたいで,要求されるままに手動で入力。ちょっと面倒だったな。今度学生に準備させるんだけど,事前に練習しておいてよかったよ・・・。

最後に Ubuntu。久しぶりに起動したので,システムそのもののアップデートをして,インストール。何だかうまくいかなかったんだけど,最終的にソースからコンパイルするとうまく走りました。Windowsより簡単ではあったな。

 

ともあれ,よいこのみんなは締め切りが迫っている中,環境を変えるのは良く無いってことを覚えておこうね!



読書会を始めます

学内で,読書会を始めます。

[amazonjs asin=”4254122144″ locale=”JP” title=”はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―”]
これをみんなで読んでいきます。

Q&Aを読み終わって,自分も何か始めないではいられなかった。半ば強引に,ゼミ生を誘い込みました。
私の考えでは,ベイズ統計学になると,もっと情報仮説が増えてくるはず。そうでなければならなかったんだ。そして,そういう仮説を思い付くには,帰無仮説を立てる訓練をしたように,統計的枠組みの中で仮説を発想する力が必要になる。

この本にあるように,発想を自由にできるようになれば,そのあとのアプローチは伝統的仮説検定の枠組みよりもずっと簡単なはず。
この読書会ではRプログラムを走らせたりすることをメインにはしない。ゆっくり,データとともに本を読んで,考える時間を持ちたいです。
一時間と短いですが,毎週ちょっとずつでも,できるだけ前に進んでいきたい。

読書会フライヤー



Lesson10.Bayesian Modeling 2

#2016.01.07現在,Rpubsが不調なため(?),こちらにアップします。

#2016.01.22 Rpubsが復旧したので,そちらにリンクを貼るようにします。なお,前に載せてあった記事にはstanコードにバグがありましたので,修正しました。

Lesson10. Bayesian Modeling 2 http://rpubs.com/kosugitti/145116



ベイズ塾例会にて発表してきました

ベイズ塾のミーティングがありました。

一人一つ,担当の「分布」を決めて発表するというお話。私はコーシー分布を担当。たまたま先日読んだ本にコーシー分布のことが書いてあったので,話が膨らんでよかったわ。

資料を公開しておきます。あんまり役に立つかどうかわからんけど。



2群の平均の差について,あるいはGitHubデビューの件

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門に打ち震えたという話から,学部生の統計の授業でもベイズを紹介したくなった。

ということで,とりあえずt検定をベイズ流にはこうやりますよ,というコードをこの本の6章を参考に書いてみた。

ただ書き写すだけでは面白くないので,ここら辺でかねてより気になっていたGitHubでの関数公開といこうじゃないか。Rstudioでバージョン管理もできるらしいしね?

ということで,GitHubやそれとRstudioの連携などに悩まされながら,なんとか書き上げてみたものがありますので,興味がある人は寄ってってねー。

 

こちらです。>> ttest2stan関数



隠れベイジアン心の叫び

[amazonjs asin=”4254122128″ locale=”JP” title=”基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門”]

この本を読むがいい!

この本を読むがいい!

これまで訳の分からぬ伝統的心理統計に悩まされてきた同胞たちよ!この本を読むがいい!
p値は仮説が正しい確率ではないと言われ(なんじゃそれ!)

差がないとは言えないがあると言い切ることはできぬと言われ(なんじゃそれ!)

同程度のサンプリングを無限に繰り返すと母数を含んでいる割合が95%だと言われ(なんじゃそれ!)

ことごとく我々文系脳の直観に反する屁理屈を積み上げ、挙げ句の果てに目の前の愛しい人が振り向いてくれないとは言えないみたいな中途半端なアドバイスしかくれなかったあの憎っくき伝統的心理統計に悩まされてきた同胞たちよ!

我々はベイジアンだ!ベイジアンになるんだ!

情報化仮説で!ベイズファクターで!主観確率まみれの人生を歩むんだ!

 

口説き落とす確率を見積もれ!(うおおおお!)

差がどれぐらいあるのかを胸を張って言い切れ!(うおおおお!)

ここからここまでの幅に答えがあると言い切れ!(うおおおお!)

 
我々はその力を手にした!

stanで!HMCで!この本で!

全国の統計ユーザーよ、ベイジアンとなって立ち上がれ!

[amazonjs asin=”4254122128″ locale=”JP” title=”基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門”]

※これはこの本の主張するところではなく、筆者が急進的ベイジアンのふりをして読んだ唄であり、良識あるみなさんは決してRTやシェアはしないでください。



MCMCでモコモコしちゃう

MCMCって大量にサンプルを発生させる技ですが,あー計算しているなー,発生させているなー,というのが目に見えてわかるようになったら楽しいと思うんです。イメージもしやすいだろうし。

ということで,Rでアニメーション画像を作る方法はないかしら,と考えたところ,animationパッケージとimage magickというのを使えば良いらしい,ということがわかりました。この2つのパッケージとアプリが用意されている前提で,次の関数を作りました。

[crayon-5bc7040c60a2d327013639/]

与えるデータは,rstanが吐き出したstanfitオブジェクトのなかでも,必要な変数のところだけextractしたやつです。
オプションとして,stepが描画するときのサンプルをいくつ飛ばしでやるか(間引き方ですが,thinだとstanの変数とごっちゃになるので名前を変えました),密度関数の線を引くかどうか(dens=TRUEで書きます。デフォルトはFALSE),スタートはバーンアウト後100回目のサンプルからにしてます。

試しに,世界一簡単なrstanコードを使って実際に使用してみたいと思います。

[crayon-5bc7040c60a40427931374/]

[crayon-5bc7040c60a45732568204/]

[crayon-5bc7040c60a50990742459/]

[crayon-5bc7040c60a55408092223/]

[crayon-5bc7040c60a5a410181568/]

最後の行でextractsしてます。

さて,このサンプルをアニメーションgifにします。
animationパッケージをrequireして,オプションをいくつか。intervalはパラパラアニメのスピードですが,0.06ぐらいが丁度いいようです。loopオプションは反復するかどうかで,デフォルトではgifが無限ループするので1にしておきます。この辺はanimationパッケージのsaveGIFオプションなので,ヘルプ見てください。

ここでsaveGIF関数の最初に,上で用意した自作関数を入れてやると,どんどんヒストグラムが育っていきます。

[crayon-5bc7040c60a60141492495/]

[crayon-5bc7040c60a64814574468/]

これを実行して得られる絵がこちら。(動かない人は画像をクリックしてね。)

MCMCexample1

モコモコしてる!MCMCがMCMCしてる!

・・・しかしここまできたら,ggplot2パッケージをつかってもう少し綺麗に描きたいな。
ということで,書き直したのが次のMCMCmovie2関数。

[crayon-5bc7040c60a6a969329965/]

使い方として,引数にchain数を入れるようにしてます。そうするとchainごとに色分けしてくれる。chainの区別なしでいいや,という人は何も書かなければ全部を1サンプルとして描いてしまいます。

[crayon-5bc7040c60a6f186542555/]

[crayon-5bc7040c60a74558524987/]

これを実行して得られる絵がこちら。

MCMCexample2

モコモコしてる!MCMCがMCMCしてる!(;´Д`)ハアハア

ちなみに棒グラフにする(dens=FALSE設定)と(なんか警告が出てうるさいですけど)

MCMCexample3

モコモコしてる!MCMCがMCMCしてる!(;´Д`)ハアハア(;´Д`)ハアハア

・・・誰かの何かのお役に立てば。



ベイズ塾で発表してきました

ベイズ塾,MCMCワークショップ第二弾で,RStanの導入について話してきました。

まあ午後の岡田先生のお話がメインなので,その前座です。

そのなかで,世界一簡単なrstanコードについて紹介したのですが,岡田先生の方から二点チェックが。

一つは,分散を逆ガンマ分布で推定するより半コーシー分布でやったほうがいい,という話。先生のお話の中で,どうして逆ガンマなのか,どうして半コーシーの方がいいのか,というのもお聞きして,納得。

あと,stanは正規分布のパラメタを平均と標準偏差で書く。普通は平均と分散でかくから,間違いガチだよね,というフォローもいただきましたが,書き直しました。ちなみにBUGSは逆数で書くのでもっと面倒だという話も。

以上二点,記事の方も修正しておきました。

R stan導入公開版 from Koji Kosugi

いろいろネタ画像を含めたのだけど(ネタ音声も),会場爆笑とならなかったのは残念です。滑り芸じゃないよ。

夜の懇親会で,次の話の機会が決まったみたいなので,次は滑らないように頑張ります?




top