Kosugitti's BLOG

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

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

stan

作戦名;Youtuberになりたい

目に見えるものはなんでもデータだと思いたい。データがあれば分析したい。

データ分析は怖くない。軽い気持ちで始めればいいんですよ。

 

 

ということが裏テーマになってます,「ポテトがあれば測りたい」。ベイズ推定の練習にどうぞ。データ,ソーススクリプトもOSFで公開中。



春の合宿セミナーで喋ってきました

行動計量学会の春の合宿セミナーは,院生の頃からずっと世話になっています。
新しい統計技術をいち早く,詳しく紹介してもらえるし,配布される資料がとても重厚。学生は参加費を低く抑えてくれているから行きやすく,至れり尽くせりという感じ。
職に就いてからも毎年情報をチェックしていて,初心に帰って勉強するつもりでよく出かけています。

が,まさかそのセミナーに講師として人前に立つ日が来るとは思わなかったなあ。
ということで,思い入れのあるセミナーで講演できました。大変うれしいことです。思い入れがあるだけにプレッシャーにも感じていて,二泊三日のブートキャンプ,スライドの枚数が400枚ぐらい準備して。
できるだけわかりやすく,教えるべき要点を絞って,何度も伝えて,とこれまでの教師経験を全てぶつけたつもりです。わかりやすかった,と受講生にいってもらえたのが何よりのご褒美でした。
受講生の皆さんが優秀すぎて,順調に進みすぎたため,スライドが足りなくなるんじゃないかと二日目の夜に少し追加したぐらい。時間調整もちょうどいい感じでできました。

埼玉の奥の方だったので移動がちょっと疲れましたが,一仕事終えて帰ってきて,心地よいバーンアウト感です。

今度はまた受講生として参加したいなあ。



Plate notation via Graphviz/RStudio

はじめに

いやー,ついに出ましたね,実戦ベイズモデリング
Amazonの書評にも投稿しちゃったけど,これプロの心理学者は全員目を通さないといけないんじゃないの。基礎実験でやってるようなことが全部ベイズモデリングされてるよ。少なくともこれからプロになる人は絶対読まないといけないんじゃないの。臨床家や実践家はなおさら。

もちろんこの本はぶっちぎってますからね,世界の最先端を垣間見せてくれるタイプのやつなので,「いくとこまで行ったらこんなこともできるのか!」とすごく動機付けられてから,青本緑本に戻って行ったらいいと思うんですけどね。

プレート表現

この本の特徴の一つが,モデルをプレート表現していること。プレート表現(plate notation)はモデルの図示にあたっての一つの表記法。怖い本で使われている表現方法で,「あー,そういうのもあるのか。そのうち勉強しないとなあ」と思っていましたが,今回この本でわかりやすく(たった2ページ!)にまとめてくれていました。これはありがたい。

これなら俺にも描けそうだ,と思ってちょっとやってみました。
というのも,以前からRStudioに標準装備のGraphvizをつかったダイアグラムの表現の仕方というのがあるのね。
詳しくはこの記事とかこの記事を見てもらえればいいかと思います。が,要するに最新版のRStudioを入れていたらもう入っているから,特に苦労しなくてももう使える環境があるってこと。あとはこのモデル図を書くdot記法とやらを習得すればいい。

RStudioが公式認定しているだけあって,ファイルの拡張子を.dotにするとdot記法の予約語はハイライトされるから,見よう見まねでいけそうなもんだ。

ということで,早速やってみました。まずは感銘を受けた付録の「プレート表現のルール」。これを次のようにして書いてみる。

plot of chunk unnamed-chunk-1

ここではDiagrammeRパッケージを読み込んでgrViz関数でdot記法のソースコードを実行させたけど,実際はこのコードを「plateRule.dot」とでも名前をつけて(拡張子を.dotにして)保存したら,RStudioのpreviewボタンを押すだけでViewerの画面に図が出て来ます。コードハイライトもされているからバッチリ!やったね!

かける!かけるぞ!

なるほどこいつぁ簡単だ,と思っていくつかやってみました。あんまり綺麗じゃないけど,セクション18.1のモデルはこんな感じ。

plot of chunk unnamed-chunk-2

もうひとつ。セクション18.11の例。

plot of chunk unnamed-chunk-3

うーん,ちょっと思ってるのと位置関係が違ったりするけど,まあなんとかなるかな?

すぐに壁が現れた

というとまあ簡単そうに見えますが,実はすぐに壁が現れたのです。上の例は,この壁の存在を理解してから描けそうなサンプルを選んだのね。

実はまず最初に書こうとしたのは,付録の図A.1なのです。一番簡単でわかりやすく表現されているから。ところがこれがうまくできない。反応の二値データは項目の特性と人の特性を分けて推定するので,それぞれにネストされているので,Graphvizでいうところのサブグラフがオーバーラップしているんですが,これが表現できない。

こうなっちゃう,

plot of chunk unnamed-chunk-4

アイテムの下にpとyijが入って来ないといけないのに,それができないんだなあ。どうにも。色々調べて見たんだけど,どうもそういう表現はできないようです(間違っていたら誰か教えて偉い人)。

Graphvizは複雑な構造とかネットワーク関係の表記にあるから,サブグラフがオーバーラップしているというのは設計思想に反するというか,考えられてないんじゃないかなと思ったり。でもプレート表現にはすごく必要なところなので(この図ではベイジアンモデリングの良さがまるで表現できてない!),ちょっと残念な感じです。

この記事ができるまで

この記事を作るにあたって参考にしたのは,次のサイト。

まあ必要そうなところを斜め読みしながら,サンプルコードを走らせながら,数時間でこの記事に描いてあるような図は描けたので,敷居は低いと言えそうですが,次のところが不満です。

  • TeX表記ができない。数式記号を「しーた」とか入れて変換するの超ダサい。下付き文字とかも書きたいし。
  • 上にも書いたけど階層構造のオーバーラップ
  • 位置の調整に細かい技が必要そう。考えなくても描いてくれるという良さでもあるんだけど。

これらの問題を解決しようとすると,今度はTeXで図を書くTikZの世界とかが出て来て,あーそこまで深入りしてらんねえ,ということで今日の自習時間は終わり。

余談だけど,もう一つの表現方法としての犬4匹本記法というのもあります。訳あって,個人的にはこっちの記法の方が慣れ親しんでる。分布の形もわかるからいいと思うんだ。

で,この記法については犬4匹本御本尊のこのサイトにあるように,TikZのソースコードの他にLibreOfficeのテンプレートがあるから,こっちのほうが導入のハードルが低いと思います。

私は当面は,ドロー系ソフトで行こうかな,と思ってます。直感的でわかりやすいからね。

ということで今日はここまで。



StanなMCMCで探索的因子分析やってみた

この記事はStan Advent Calender 2016にエントリーしたものです。

やっていることは,ここ数年関わっている探索的因子分析をベイズ推定するっていう話なのですが,最初に言っておきますと,探索的因子分析ってベイズモデリングに向いてないと思うのですよw

ベイズモデリングが何か,という話は一言では終わらないですが,私は「データ生成メカニズムをデザインするモデリング」と理解しています。なので,認知心理学的なモデルなんかにはもってこいだと。そもそも心理学は個人差をどう考えるか,というような世界でもあるので,階層的なモデリングのためにとてもいいと思うのですね。

一方,探索的因子分析はその名の通り「探索的に」データを捉えるわけですから,データがどういう構造しているかわからない,というのはそもそもちょっと違う。もちろん,n因子構造だと仮定してモデルを書くと,ということになるんですが,CFAならともかくEFAだとまあありとあらゆる係数を推定していくことになるのです。なので,ちょっと大変だなあ,と書き終わってから思うたわけです。

 

さてさて,本題です。
因子分析というのは,$$N$$人の個人$$i$$が$$M$$個ある項目$$j$$に反応したデータセット$$x_{ij}$$(ただし$$N \ll M$$,つまり項目数より個人の数が十分に多い)について,$$k$$個の因子という潜在変数を仮定して,因子と項目との関係を表す因子負荷量$$a_{jk}$$と,因子と個人との関係を表す因子得点$$f_{ik}$$との積和の形で表すモデルです。

$$ x_{ij} = \mu_j + \sum_{k=1} a_{jk}f_{ik} + e_{ij}$$

ここの$$e_{ij}$$は誤差項ね。これをベイズ的に考えるときは,基本的に潜在変数が正規分布に従うことを仮定する,というだけなので,それをそのまま素直にモデリングしてあげれば一丁あがりです。

ただ,探索的にやる場合は因子負荷行列に制約をかけてあげる必要があります。因子負荷行列は一般に,直交行列になるような制約を置いてあげるものですが,SEM的には$$a_{jk}=0$$(ただし j < k の場合)という制約をかけることになります(豊田先生の「共分散構造分析【技術編】」を参考)。これをコード化したら次のようになります。

ここのtransformed parametersブロックで面倒なことをしています。整数でindexというのを宣言して置いて,このループの中だけでカウントするインデックスを用意しました。もう少しうまいやり方もあると思うんだけど,とりあえずこれでできます。

 

ただし!このやり方はすごく時間がかかります。時間がかかる理由は,データサイズに比例して推定しなければならないパラメタが多いということ。N人A因子あればN*A分の因子得点を推定しますからね!こりゃめんどうです。

実際,一般的な因子分析は,基本モデルにいくつかの条件(因子得点は標準化得点,誤差の平均はゼロ,誤差因子と共通因子の相関はゼロ)を加えて,

$$ \Sigma = A \Phi A’ + \Delta $$

特に直交因子モデルの場合は$$\Phi$$は単位行列なので

$$ \Sigma = AA’ + \Delta $$

とした分散共分散行列で考えることになります。

Stanで多変量正規分布を推定するのはなぜか大変重たいので,ウィシャート分布から推定することにし,因子負荷行列もこれスキー分解したベクトルで与えてやると,Stanコードは次のように大変短くなります。

データから平均値を取り去った(中心化した)データセットとして与え,それの積和行列がtransformed dataブロックでspdとして作られており,それがwishart分布に従うようになっています。

これでStanはちゃんと推定してくれます。が!これでお話は終わりではありません。実はここから得られるMCMCサンプルの推定値を考えても,0ばっかり!みたいなことになるのです。それには三つほど問題が関係しています。

一つはMCMCサンプルごとに出てくる因子の順番が違う問題。これはサンプリング後に因子の順番をソートして整えることで対応できます。

もう一つはMCMCサンプルごとに出てくる因子負荷量の符号が違うこと。これは因子毎の負荷量の総和が正になるように整えるとか,任意のMCMCサンプルをつかってその因子毎の負荷量の総和の符号に合わせる,といった恣意的な基準を作って対応できます。

もう一つは因子の回転についてです。stanの中では回転作業をしていませんので,得られたMCMCサンプルを全部回転することで対応できます。

 

・・・とここまで書いておわかりのように,サンプルを得た後の作業が結構大変なんですね。実際にやってみると,サンプリングよりその後の回転やソートのほうが時間がかかったりします。

その作業をするコードを載せても良いのですが,それはRの話でStanの話じゃなくなってしまいますから,この記事はここまでです。Rの関数として,psychパッケージのようにfa,と書いたら探索的ベイジアン因子分析してくれるものを現在開発中ですので,興味がある方は個別にご連絡を。

 

それにしても,確率モデルをほぼそのまま書くだけでできちゃうStanって,本当に便利ですよね。みなさんもEnjoy Stan!

参考文献;

[amazonjs asin=”4254126646″ locale=”JP” title=”共分散構造分析 技術編―構造方程式モデリング (統計ライブラリー)”]



習作;ベイズファクター

岡田(2014)のベイズファクターの計算をStanでやってみた件。

これについては,次のスライドに詳しい。

ただ,途中のStanコードがちょっとわかりにくいのと,計算をstan内部でgenerated quantitiesを使ってやったらいいんじゃないかと思ったので,そのように書き直して見た習作。

推定結果は論文とほぼ一致した。BF1,BF2,PMP1,PMP2を参照。



Stanで体重の推移をみつめてみた(状態空間モデル)

この記事は Stan Advent Calendar 2016 のエントリーです。

Stanいいですね。いいですね,Stan.

私はもっぱらRから触るrstanユーザーなんですが,いろいろなところから触れる間口の広さがいいですね。

もちろん「尤度関数を書かなくても推定してくれる」とか「細かい設定をしなくても大体ちゃんとやってくれる」といったStan本来のありがたみも忘れちゃいけません。

なにより,ベイジアンモデリングのハードルがすごく下がった気がします。頭の中でイメージしていることが,そのまま書けるっていうのは,とてもいいことだと思うのです。

この「頭の中のイメージがそのまま書ける」ということを一番実感したのは,状態空間モデルのお勉強をしているときでした。なので,その時のことを記事にして,Stanへの感謝とエールに代えたいと思います。

時系列データと状態空間モデル

時系列的なデータの場合,一つ前の時点との関係(自己相関)が強いもんだから,そのままの分析じゃダメで,時系列分析という別の考え方が必要になる。これはまあいいです。

でも,従来の(最小自乗法,最尤推定法的な)時系列分析はなんだかわかりにくい。

時系列分析については,いろいろな本も出ています。あるいはこちらのサイトはすごく丁寧にわかりやすく書いてくれていますので,時間かけて読めば何をやっているかはちゃんとわかります。でも,ちゃんと面倒なことがわかります。

というのも,とある時点での推定値がでても,それが期間全体で成立するかどうかを前後にわたってチマチマ調整していかなければならないからです。私の理解したイメージでは,基本的にt時点からt+1時点への回帰で,それがワンステップ前,ワンステップ後でも推定地として大丈夫かどうか,均してからつかうという感じ。すごい力技な感じがするんですね。実際,モデルを決める,パラメタを求めて,フィルタをかけて,スムージング,という順番でやらないといけないのです。

そこで出ました状態空間モデル。これは

  • 目に見える観測の世界
  • 目に見えない状態の世界

を区別して,状態は連綿と続いていく,観測はそれに誤差が乗っかる,というイメージでモデリングしていくのですが(これについてもこちらのサイトが詳しいです!),Stanはこのイメージをモデルに書くときに,本当にそのまんま!という感じなのです。これを実際の例で見ていきましょう。

 

状態空間モデルによる体重の予測

私はなぜか,毎朝起きたら体重と体脂肪を測定して記録する,というのを日課にするようになっているのですが,この習慣が2012年ごろから続いており,時折記録できない日(出張で自宅にいない=いつもの体重計がない場合)なんかもありつつ,2016年7月末の段階で1077レコードありました。

実は昨年末帰省したときに,京都ラーメンが美味しすぎてたくさん食べていたら,体重が最高83kgまで達しまして,これはまずいんじゃないか(当方身長173cm)と,今年に入ってから少し体重を下げる努力をしております。この辺の効果なんかも考えてみたいところです。

 

さて,状態空間モデルは「観測値」と「状態」を区別すると言いました。ここでは,観測値=体重計の目盛りです。でもうちの体重計は,結構数字を出すのに悩むんですね。77.4?77.5?77.3?いやまて77.7・・・?とか考えてから,今日は77.6kgです,と表示したり。体の傾け方とかで数値がぶれるのかしら,と思うのですが,まあこの目盛りは信じるしかない。

こうした測定誤差とは別に,肉がついてるとか落ちたとかいう実際の状態もあるわけで,そこを別なものとしてモデリングします。下のイメージ図参照。

状態空間モデルのイメージ図

さて,これをStanでかくと,モデルの部分はまさにこのイメージのまま,なのです。

いかがでしょうか。観測データ(W[n])は実際の肉(mu[i])に誤差(sdV)がついて変化します。実際の肉(mu[i])は昨日の肉(mu[i-1])に誤差(sdW)がついた形で変化します。
モデルのイメージ図をそのまま書いただけになっているのです!フィルタとかスムージングとかいりません!

欠測と予測

自宅の体重計で毎朝測るのですが,自宅で目覚めない日もあるわけで,そういう日は測定できません。欠損値になります。

状態空間モデル2欠測

それでも肉がなくなるわけではないので,状態は続いていきます。これをモデリングする場合は,測れなかった日の観測値もパラメータとして推定してやれば良いのです。

Stanは残念ながら欠測値に対応していませんから,「これは欠損だよ」「観測データはちゃんと揃ってるよ」という形で渡してあげなければなりません。私はここで,欠測のところに特殊な数字(9999)を入れて,Stanに渡してあげることにしました。

ついでに。
時系列的な変化は,やはり未来に向けての予測がしたいですよね。予測,というのも「(まだ)観測されていないデータ」という意味で,欠損値ですから,同じアルゴリズムが使えます。
状態空間モデル3予測

予測したい日数を欠損としてデータセットに加えてやれば,この問題は解決です。先ほどのコードにちょいと付け足します。

これでよし。さて,Stanのコードはどうなっているかというと,

観測方程式のところに,missを埋めていくインデックス,nmissというのを作りました。これは0からカウントしていきます。

もし観測(w[i])が9999じゃなかったら,それまでどおり肉の状態(mu[i])に誤差がついた形であらわれているのですが,もし9999だったら=欠測だったら,肉の状態(mu[i])からでてきた推測すべきパラメータ(Miss_W)になるよ,ということです。

そして欠測していようがしていまいが,肉の状態は1-step前の状態から推移し続けます。

これで推定してみた結果は次の通り。

 

上が観測の誤差成分,下が状態の誤差成分です。観測成分による誤差が状態の誤差の約2倍ぐらいあります。
実際のお肉の状態は次のような感じ。

ほうほう・・・。

さて,この記事は事前に描いて今日公開するよう予約投稿してあるのですが,今年のお正月から7月末までの前半期をつかって,この12月1日の体重を予測しますと!

実は74.1950kgなお肉の状態になっているはずです!当たってるかい?今日の俺よ。
まあ95%確信区間でいくと69.86-78.03kgの間に入って入ればStanの勝ちということです!
確信区間も含めて推移をグラフにしてみると,測定が終わった次の日からパッと幅が広がっていくのがよく見えます。そりゃそうだよねー。

 時系列予測

ともあれ,いわゆる時系列分析をするよりも,Stanでもベイジアンモデリングしたほうが直感的にわかりやすいし,簡単だと思うので,超オススメです。
もちろんこの「状態」のところに,予測するモデルを入れていって,より複雑なモデリングにすることができます。私は曜日によって違うんじゃないかとか,出張の前後で体重が変わるんじゃないかとか,いろいろ細かな変数を入れて推測させて見ましたが,ほとんど効果がなかったので,ここではそうしたモデリングの部分はご紹介しませんでしたが。

 

 

さあEnjoy Stan!



不倫報道の間隔を推定する

2016年は不倫報道が多いですね!全くどうでもいいんですが,これほど多いとデータとしてみなし,なんかの分析をしたくなる。

そこでデータをつくって分析してみます。

「報道が生じる」のがポアソン分布,ポアソン分布で次の観測が生じるまでの感覚は指数分布に従うと考えられるので,これをモデル化。

結果のmuの値が示すのが,次の不倫報道まで待つべき日数になります。

もちろんこれは,不倫「報道」までの日数予測であって,かなり人為的な要素があることです。予測の精度を検証しようってな話ではありませんが,占いによる予想よりはよっぽどマシでしょう,という笑い話としてお納めください。



Hijiyama.R #5でプレゼンしてきました

こんかいのHijiyama.RはRをつかって文書作成,R markdown祭りという裏テーマがありまして,それに乗っかって話をしてきました。

内容は,心理統計の課題レポートを課すとき,乱数でデータ作るんですけど,思い通りのものができませんよね・・・というところから始まって,実はベイジアンになろうね!という話になっています。うふふ。

2016.10.17追記

いくつかソースをアップしておきます。

今回のデモで使ったRmdファイルはこちら。>> ReportExample.rmd

出来上がるPDFファイルはこちら。>>ReportExample.pdf

後半出てきたStanコードはGistに上げました。これらのコードはあくまでも上の課題に対応した表記になっている(ハイレベルスクリプトではない)ので,自分の環境に合わせて書き直したり,別途キックするRコードを書いたりする必要があります。

対応のないt検定のstanコード

対応のあるt検定のstanコード

一要因群間ANOVAのstanコード

一要因群内ANOVAのstanコード

混合計画ANOVA(間2,内2)のstanコード

回帰分析のStanコード



習作;ラベルスイッチング問題をなんとかする関数

前回の続き。Stanコードに対数尤度を出す部分を追加。。。
(どう見てもカッコ悪い。スマートに書き直したいものだ。ていうか,そもそもこれで合ってるのかな?)


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{
  real log_lik[N];
  simplex[K] u[N]; #class membership probability

  for(n in 1:N){
    u[n] = softmax(ps[n]);
    log_lik[n] = log_sum_exp(ps[n]);
  }

}

で,鎖が複数になるとラベルスイッチング問題が生じるんだけど,これを解決するlabel.switching関数を使ってみる。

## どこに分類されたか
allocK < - rstan::extract(fit_sp,pars="u")$u
# 各チェインごとに割り当てられたクラス
zChain <- matrix(ncol=N,nrow=itr*C)
for(i in 1:(itr*C)){
  zChain[i,] <- apply(allocK[i,,],1,which.max)
}

## MCMCで得られたパラメタ
mcmc.params <- rstan::extract(fit_sp,pars=c("mu","sig2","theta"))

# mcmc * num.of.clusters * num.of.params
mcmc.pars <- array(data=NA,dim=c(itr*C,K,3))
mcmc.pars[,,1]<- mcmc.params$mu
mcmc.pars[,,2]<- mcmc.params$sig2
mcmc.pars[,,3]<- mcmc.params$theta

# この関数はSJW法を使う時に必要になってくる。completeオプションに関数を渡さないといけないので!
complete.normal.loglikelihood<-function(x,z,pars){
  #	x: data (size = n)
  #	z: allocation vector (size = n)
  #	pars: K\times J vector of normal mixture parameters:
  #		pars[k,1] = mean of the k-normal component
  #		pars[k,2] = variance of the k-normal component
  #		pars[k,3] = weight of the k-normal component
  #			k = 1,...,K
  g <- dim(pars)[1] #K (number of mixture components)
  n <- length(x)	#this denotes the sample size
  logl<- rep(0, n)	
  logpi <- log(pars[,3])
  mean <- pars[,1]
  sigma <- sqrt(pars[,2])
  logl<-logpi[z] + dnorm(x,mean = mean[z],sd = sigma[z],log = TRUE)
  return(sum(logl))
}


# パッケージのお出まし
library(label.switching)
# 全方法試す
set <- c("STEPHENS", "PRA", "ECR", "ECR-ITERATIVE-1", "ECR-ITERATIVE-2", "AIC", "SJW","DATA-BASED")

# 対数尤度
log_liks <- rstan::extract(fit_sp,pars="log_lik")$log_lik
library(loo)
loo(log_liks)
waic(log_liks)

#ピボットは対数尤度が一番大きいやつにしてみた(なんでもいいと思う)
# MCMCサンプルごとに対数尤度が一番大きいやつを算出(apply部)
# 表の形で集計,大きい順に並べかえる(sort(table)部分)
# 表のラベルを取り出して(names部分),数字にする(as.numeric)
pivot <- as.numeric(names(sort(table(apply(log_liks,1,which.max)),decreasing=T)[1]))
# さあ実行
ls <- label.switching(method = set,
                      zpivot = zChain[pivot,],
                      prapivot = mcmc.pars[pivot, , ],                      
                      z = zChain,
                      K = 3, 
                      p = allocK,
                      mcmc=mcmc.pars,
                      complete=complete.normal.loglikelihood,
                      data = X)
# 各推定法で結果が一致したかどうか
ls$similarity
# 推定された所属クラスタ番号
ls$clusters

練習ここまで。




top