The One with ...

研究用のメモなどを公開しています

ggplot2による事後分布の視覚化

絶対忘れるので自分用にメモしておく.

目的

Stanで推定した結果を分布のまま直接比較したい.

ggplot2用のデータフレームの作り方

流れはこんな感じ.

  • Stanで推定した結果をfitに記録(CSVで記録すると結果的に遅かった.sample_file からのサンプル列の取り出しは面倒.たいした量じゃないので今回はRDataを保存してloadする).
  • stanfitオブジェクトをas.data.frame()でデータフレームに変換
  • パラメータの各列を1列にスタック(95年と2015年で年齢幅が違うのでindexを小細工したうえでforを回した)
  • ggplot用のインデックス列を作成(年齢,性別,調査年).インデックス列は同じラベルがmcmcサンプルの数だけならぶ.今回は1000×3(chain)個.
  • このときスタックの順番は男女男女とし,性別のインデックスを行方向に対応させる.
  • plotしたいパラメータとインデックス列をデータフレームにまとめ,ggplotにわたす.
  • ggplotのfill=で1つの枠内に男女を並べ,facet_wrap(~ age, ncol = 5)で年齢別に並べる.結果年齢別に男女比較が行列型に並ぶ
  • 調査年比較は,パラメータをスタックする際に,2015,1995の順番で並べる.

コードはこんな感じ. fitに2015年度データの推定結果,fit2に1995年度データの推定結果が記録されている.

結果の一例.

f:id:hamada7418:20170817150851p:plain

暇があったら,tidyverseでもうちょっと可読性を高めたい.

多次元正規分布のパラメータ推定

今日はStanを使って多次元正規分布のパラメータを推定する練習をやってみましょう(授業の続き).

まず推定用のデータをRで生成します.

library(MASS) # for mvrnorm( )
N0 <- 300 #データ個数
m1 <- 0;m2 <- 0;s1 <- 1;s2 <- 1;r <- 0.6 #パラメータ設定
mu = c(m1,m2)  #平均ベクトルの定義
sigma = matrix(c(s1^2, r*s1*s2, r*s1*s2, s2^2), 2, 2)#共分散行列の定義
#r=(r*s1*s2)/(s1*s2) 共分散行列は基準化しているので,rは相関
data0 = mvrnorm(N0, mu, sigma)
#2次元正規乱数をN0個を生成  相関cor(exv[,1],exv[,2]) 

2次元の正規分布にしたがう乱数を発生させる関数はmvrnorm(N0, mu, sigma) です(MASSパッケージの関数).引数のN0が個数,mu, sigmaは平均ベクトルと共分散行列です. 相関rをダイレクトに指定するため基準化しています.

この平均ベクトルと共分散行列がStanで推定すべきパラメータです.正しく推定できれば,自分が定義した 真のパラメータに近い値が出てくるはずです.

次の図はmvrnorm( )を使って生成したデータの散布図です.相関が0.6なので,右方向に点が上がっていく様子が見て取れます.多次元正規分布を使えば, このように指定した相関を持つデータ列を簡単に作ることができます.

f:id:hamada7418:20170712163246p:plain

コード

R上の分析実行コードですStanコードは例によってアヒル本記載のコードを流用しています.

Stanコードの中で注目すべきは,共分散行列用の型であるcov_matrixという型を使っている点です. 具体的には,Stanのパラメータブロックが次のように書かれています.

parameters {
  vector[D] mn;
  cov_matrix[D] cov;
}

この型は,対称な正定値行列に対応しています.

Stanにはそのほかにも単体用のsimplexや, 順序を持ったベクトルに対するorderedなどの型も用意されています.

さて推定結果fitの中身を確認してみましょう.

            mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
mn[1]       0.03    0.00 0.06   -0.09   -0.01    0.03    0.07    0.14  2853    1
mn[2]      -0.03    0.00 0.06   -0.14   -0.07   -0.03    0.01    0.08  2660    1
cov[1,1]    1.05    0.00 0.08    0.90    0.99    1.04    1.10    1.22  2852    1
cov[1,2]    0.59    0.00 0.07    0.47    0.55    0.59    0.63    0.73  2640    1
cov[2,1]    0.59    0.00 0.07    0.47    0.55    0.59    0.63    0.73  2640    1
cov[2,2]    0.98    0.00 0.08    0.84    0.93    0.98    1.03    1.14  3000    1
lp__     -236.19    0.04 1.53 -239.99 -236.92 -235.88 -235.09 -234.19  1367    1

Samples were drawn using NUTS(diag_e) at Wed Jul 12 16:05:24 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

おおむね,うまく推定できています*1

plot(fit)で95%ベイズ信頼区間を確認しましょう.

f:id:hamada7418:20170712173334p:plain

stan_hist(fit)でパラメータのヒストグラムを描きます.

f:id:hamada7418:20170712173337p:plain

パラメータの真値( m1 <- 0;m2 <- 0;s1 <- 1;s2 <- 1;r <- 0.6 の部分)をいろいろ変えながら 自分で試してみましょう.

解析的な事後分布の導出

さてここからは事前共役分布を使った事後分布の導出について話します.まず1次元の場合,{ \displaystyle  \mu   }に対する事前分布を{ \displaystyle
\mu \sim N( \mu_{0}, \sigma_{0}^{2}) }とすると,事後分布のカーネル

{  \displaystyle p ( \mu | y ) \propto \text{exp} \left\{ - \frac{\mu - \mu^{*})}{2 \sigma^{*2}} \right\} }

となる.ここで

{  \displaystyle \mu^{*}=w \mu_{0}+(1-w) \bar{y}, \sigma^{*2}=\frac{1}{1/\sigma_{0}^{2} + n/ \sigma^{2}} } である.

またウェイトは

{ \displaystyle w = \frac{1 / \sigma_{0}^{2} }{1 / \sigma_{0}^{2} + n/ \sigma^{2}} }

である.つまり \mu の事後分布は正規分布となる.

(続く.はてぶろのTeX表記,エスケープが多すぎて安定しない…)

*1:当たり前やん,と思うかもしれませんが,「正解」が存在する推定の問題を自分で作っているのだと 考えてください.もちろん経験的には意味はありません.しかし分析手法を理解するために,これは最適な方法の一つです

回答のランダム化とStanによる離散パラメータ推定

Stanによる離散パラメータの推定

離散パラメータが含まれている場合のMCMC推定は基本JAGSでやってきたのだが,新しい方法を勉強するためにStanで周辺化消去を試してみる. 以下StanとRコードは,松浦さんのアヒル本11章に依拠しています.

回答のランダム化

被験者が答えにくい質問への回答の母比率を推定するための手法に,「回答のランダム化」という方法がある.これは次のような方法だ.

まず被験者がコインを振る.調査者はその結果を見ない.コインの裏が出たら,被験者は質問に常にYESと答え,表が出たら真実の回答を述べる.自分が真実の回答を述べたか,コインによってYESと答えたのかを,調査者に知られることはないので,被験者は安心して回答できる.よく工夫された設計である.

この調査法は規範的・倫理的に真実の回答を表明しにくい質問に適している.例えばローゼンタールのテキスト『運は数学にまかせなさい』では,マリワナ喫煙者の割合の推定例としてこの手法を紹介している.

人工データの作成

推定がうまくいったかどうかは,そもそも真の分布や真のデータ生成メカニズムが特定できていないと判定のしようがない.そこでまずR上で推定すべき「真のデータ」を生成する.

Rにはベルヌーイ分布の確率変数を生成する関数がどういうわけか用意されていないので,まずは関数bern( )を作成し,それを使ってデータを作る.

人工データを作成する部分の説明

まずbern(0.5)を使い,0 or 1の乱数を発生させる.0がでたら必ず1を出力し,1が出たら確率q0で1を出力(1-q0で0を出力)する.これを1000人分繰り返してデータYとして保存する.

Stanコードはアヒル本記載のコードをそのまま使っています.この例では推定すべき真のパラメータ,つまり母比率はq0=0.85であると仮定している.

Stanの結果は以下の通り.

> fit
Inference for Stan model: model11-1.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
q       0.86    0.00 0.02    0.83    0.85    0.86    0.87    0.89  2636    1
lp__ -253.67    0.02 0.73 -255.81 -253.80 -253.40 -253.22 -253.17  1439    1

Samples were drawn using NUTS(diag_e) at Fri Jul 07 16:48:00 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

結果qの中央値と平均が0.86で,95%ベイズ信頼区間は(0.83~0.89)となった.

確率論による母比率の(解析的)推定

1000人の回答者に「マリワナ喫煙の経験がありますか」と質問して,コインによってランダム化された回答データがあるとしよう.表が出た場合のみ本当の回答をして,裏が出た場合は常に「はい」と答えてもらう.このとき,MCMCで推定しなくても,母比率の点推定は可能である.

考え方はこうだ.まず集合を4つ定義する.

 A=\{ コインで表が出た人\}

 B=\{ コインで裏が出た人\}

 C=\{ Aの中で実際にマリワナを吸った人\}

 D=\{ Bの中で実際にマリワナを吸った人\}

データを集計したところ,質問に「はい」と回答した人数が930だったとする.つまり,

|B \cup C|=930

である.また大数の法則により,コインは大体半分は裏が出ていると考えられるので|B|=500と仮定する.同様に |A|=500と仮定する.またコインの裏表がランダムに出ているから,

 \displaystyle \frac{|C|}{|A|}=\frac{|D|}{|B|}

と仮定できる.要するに,表が出た人々における「実際のマリワナ喫煙者割合」と,裏が出た人々における「実際のマリワナ喫煙者割合」は等しいはずである.

この値は,全体における「マリワナ喫煙者の割合」に一致する.

 C の要素数が分かれば,\frac{|C|}{|A|}がすぐ分かるので, |C|を計算する.

まず

|C \cup B|-|B| =|C| \quad (B \cap C=\emptyset より)

である.つまり930-500 = |C| だから 430 =  |C|である.

よって, \displaystyle \frac{|C|}{|A|}=\frac{430}{500}=0.86
である.

この数値はベイズ推定による事後平均およびMAP推定値と一致している.

Stanによる推定例:ベルヌーイ分布のパラメータ

連勝記録から勝率を推定する

ある棋士が29連勝した,というデータから,勝利確率をMCMC推定してみる.統計モデルの仮定は次の通りである.

  • 毎回の勝ち負けは独立にベルヌーイ分布に従っている

  • ベルヌーイ分布のパラメータpの事前分布は範囲[0,1]の一様分布とする

StanとRのコードは以下の通り

パラメータpの事後分布のプロット.

f:id:hamada7418:20170628100158p:plain

ベルヌーイ分布にしたがう確率変数の例を日常生活の中から見つけてみよう.

ちなみに勝率0.9で29連勝する確率は

0.929=0.04710129

である.

解析的な事後分布の導出

上述の計算機による推定と同様の結果を,手計算によって導出してみよう.

試行数1のベルヌーイ分布の確率関数は

{ \displaystyle
p ( Y | \theta) = \theta^{Y} (1-\theta)^{1-Y}
}

だから,n回試行の同時確率は

{ \displaystyle
\prod_{i=1}^{n} \theta^{Y_{i}} (1-\theta)^{1-Y_{i}}
}

である.この同時関数をパラメータθの関数として見たものが尤度関数{p(y|\theta) }である. θの事前分布としてベータ分布

{ \displaystyle
B(a,b)
}

を仮定する.さて事後分布は

{ \displaystyle
p(\theta|y) \propto p(\theta)p(y|\theta)
}

{ \displaystyle \quad \propto \theta^{a-1} (1-\theta)^{b-1} \prod_{i=1}^{n} \theta^{Y_{i}} (1-\theta)^{1-Y_{i}}
 }

であった.ここで積の前半部分はベータ分布のカーネル 

{ \displaystyle \theta^{a-1} (1-\theta)^{b-1} } 

であり,後半部分はベルヌーイ分布の尤度

{ \displaystyle \prod_{i=1}^{n} \theta^{Y_{i}} (1-\theta)^{1-Y_{i}} }

である.その計算結果は

{ \displaystyle
\theta^{\alpha-1} (1-\theta)^{\beta-1}  }

ただし { \displaystyle \alpha = a + \sum_{i=1}^{n}y_{i}, \quad \beta = b + n - \sum_{i=1}^{n}y_{i} } 

となる.これはベータ分布のカーネルである.つまり事前分布としてベータ分布を指定すると,ベルヌーイ分布のパラメータの事後分布はベータ分布

{ \displaystyle
\theta | y \sim B( a + \sum_{i=1}^{n}y_{i},  b + n - \sum_{i=1}^{n}y_{i})
 }

となる.よって事前分布としてa=1, b=1と仮定すればθの事後分布は

{ \displaystyle
B(1+29, 1+29-29) = B(30, 1)
 }

である.このベータ分布をθの事後確率分布としてplotとしてみよう.まず事前分布がB(1,1)の場合.

f:id:hamada7418:20170713000704p:plain

次に事前分布がB(2, 2)の場合.

f:id:hamada7418:20170713000707p:plain

Stanで推定した結果に対応した確率密度関数のplotを得た*1

*1:厳密に言うと,Stanの結果のplotはMCMCサンプルのヒストグラムなので,確率密度関数のplotではない.

StanからWAICを計算するためのメモ

StanからWAICを計算するためのメモ

豊田本(2017)を参考に,Stanの結果からWAICを計算する方法のメモ.

テキスト198ページ記載のコードは略記のため,そのまま書き写して実行してもStanが動きません*1. 公開されたコードを参考に確率分布のパラメータを適切なブロックで定義すれば,Stanが動きます.

まずStanコード内でモデルの対数尤度を生成量(generated quantities)として定義します.

次にR上でlooパッケージを使って,Stanのサンプリング結果から対数尤度を取り出し,waic関数にわたします*2

それぞれのコードはこんな感じです(ポイントだけ省略して書いています.実際にはデータブロックやパラメータブロックを書かなくてはいけません).

上記のGitHub Gistからの埋め込みは,Gist上で「Embed」をクリップボードにコピーしたものを ブログのソースに張り付けるだけで完了します.

*1:モデルブロック内で宣言したパラメータが存在しないというエラーメッセージがでます.JAGSならこれでも行ける気がしますが

*2:内部的に対数尤度の総和はtargetに自動記録されるけど,平均計算用に対数尤度をベクトルで記録する必要がある,という理解でいいのだろうか?

TeXstudioでオリジナルコマンドのショートカットキー設定方法

TeXstudio便利ですね.私は長い間 秀丸+マクロ を使っていましたが,TeX Live導入後は TeXwokrsTeXstudioを併用しています.

さて,TeXstudioでオリジナルのショートカットキーを作成したいと思うことが時々あります.

時間がたつと,すぐにやり方を忘れてしまうので,未来の自分のためにメモを残しておきます.

例.ショートカットキーを使って,ワンストローク

 \begin{align*} 
 &   \\ 
 & 
 \end{align*} 

といった長い数式(環境)を入力したい,という状況を考えます.

方法の概略

「メニュー」にまず登録して,その後ショートカットキーを設定する.

手順
  • メニューから「オプション」→「TeXstudioの設定」を選ぶ.

  • 「高度なオプション」にチェックを入れる(チェックしないと「メニュー項目」が出てこない)

  • 左の列から「メニュー」を選ぶ.次に「数式」の項目から適当な一つを選び,右クリックする.

f:id:hamada7418:20170613155650p:plain

  • 右クリックメニューから「(前に)新規メニューアイテムを挿入」を選択する.新しい項目が挿入されたら,適当な名前を入力する(今回の例では「複数行」という名前にしたと仮定します).

  • 「コマンド」の列に,挿入したいコマンドを書く.例えば

\begin{align}%\ & \ %\ & %\ \end{align}

のように書くと,実行したときにソース上で

 \begin{align*}   
 &   \\   
 &     
 \end{align*} 

が挿入される.メニュー設定のコマンドでは「 %\ 」が改行を意味しています. 「 %| 」で挿入後のカーソル位置を指定できます(必要がなければ指定しなくてもOKです).

  • メニュー登録が完了したら,左の列で「キーボードショートカット」を選択.「数式」の項目を展開すると,今登録した「複数行」があるので「現在のショートカット」のセルをダブルクリックする.

f:id:hamada7418:20170613155654p:plain

  • そこで自分の好きなキーを登録(例として「Ctrl + E 」を登録したと仮定します).登録後,OKボタンを押して設定を終了する.

  • ソース上で「Ctrl + E 」を押す.

ソースに

 \begin{align*}   
 &   \\   
 &     
 \end{align*} 

が挿入されれば成功です.

うまくいったら,同様のやり方でよく使うコマンドを登録して,入力の手間を省きましょう(もともと登録されているものは,ショートカットキーを編集するだけでOKです.なおメニューに登録していないコマンドをショートカットキー登録する方法は知りません).

tikz-decoration

TikZdecorationがうまくいかないので検索してみると,次のような解決策がweb上にあった.

マニュアルに記載の

\draw[decorate,decoration={coil,segment length=4pt}] (0,1) -- (3,1);
\draw[decorate,decoration={coil,aspect=0}] (0,.5) -- (3,.5);

を実行するには,

\usetikzlibrary{decorations}

ではダメで,

\usetikzlibrary{decorations.pathmorphing}

まで指定しないといけないらしい.

マニュアル原典が大部すぎるため,これは確かにわかりにくい.