diff --git a/example1.R b/example1.R index 9c588ab..cfbcc75 100644 --- a/example1.R +++ b/example1.R @@ -57,7 +57,7 @@ post2 <- coda.samples(model2, variable.names = "lambda", ## Plot trace #pdf("example1-2.pdf", width = 360/72, height = 240/72, # family = "Helvetica", pointsize = 10) -traceplot(post2, las = 1) +traceplot(post2, las = 1, col = c(1, 2, 4)) #dev.off() diff --git a/example2.stan b/example2.stan index ec67029..7d61011 100644 --- a/example2.stan +++ b/example2.stan @@ -11,9 +11,7 @@ parameters { } transformed parameters { - vector[N] logit_p; - - logit_p = beta + beta_x * X; + vector[N] logit_p = beta + beta_x * X; } model { diff --git a/mcmc_text/example1-2.pdf b/mcmc_text/example1-2.pdf index b7c3128..2020706 100644 Binary files a/mcmc_text/example1-2.pdf and b/mcmc_text/example1-2.pdf differ diff --git a/mcmc_text/mcmc_text.tex b/mcmc_text/mcmc_text.tex index c41cb19..2c1dc58 100644 --- a/mcmc_text/mcmc_text.tex +++ b/mcmc_text/mcmc_text.tex @@ -55,8 +55,8 @@ \section{はじめに} \textsf{JAGS}\cite{JAGS}を使用する。\footnote{本テキスト執筆時には，R 3.4.1, JAGS 4.3.0, Stan 2.16.0を使用した。} %（実習には\textsf{MCMCpack}を使用する）。 また，\textsf{Stan}\cite{Stan}の紹介もおこなう。 -例題等のファイルは\url{https://github.com/ito4303/naro_toukei}で -公開しているので，参照されたい。 +例題等のファイルの最新版は\url{https://github.com/ito4303/naro_toukei}で +公開しているので，必要に応じて参照されたい。 \subsection{統計モデリング} @@ -67,7 +67,7 @@ \subsection{統計モデリング} 統計モデリングにおいては，複雑なモデルが必ずしもよいというわけではない。 現実世界の本質的な部分のみをうまくモデリングすることにより， -現実をよく理解できたり，うまく予測できたりするのである。 +現実をよく理解できたり，うまく予測できたりするのである\cite{BPA}。 \subsection{ベイズ統計モデリングとMCMC} @@ -167,8 +167,10 @@ \subsection{MCMCのアルゴリズム} \section{MCMCのためのソフトウェア} CやFortranなどでMCMCのプログラムを作成することも可能だが， 専用のソフトウェアを使う方が簡単であろう。 +MCMCのソフトウェアには以下のようなものがある。 + \begin{itemize} -\item \textsf{R}上で: MCMCpackなど +\item \textsf{R}のパッケージ: MCMCpackなど \item 専用ソフトウェア: \textsf{WinBUGS}, \textsf{OpenBUGS}, \textsf{JAGS}, \textsf{Stan}など \begin{itemize} \item 上に挙げたもののうち，\textsf{Stan}以外は，BUGS (\textbf{B}ayesian inference \textbf{U}sing @@ -216,50 +218,50 @@ \subsection{WinBUGS} で \textsf{R}と連携可能なので，そちらから使う方が使いやすい（と思う）。 \end{itemize} -\subsubsection*{インストール} - -ウェブサイトからインストーラーをダウンロード -できる。いずれの環境でも\textsf{WinBUGS}のインストール後には1.4.3パッチをあ -て，``key''をインストールすること。``key''をインストールすることにより， -全機能が使えるようになる。 - -\paragraph{Windowsへのインストール} - -インストーラー（\texttt{WinBUGS14.exe}）を起動して，あとは指示に従ってい -けばインストールされる。 -ただし，Windows Vistaでは，\texttt{C:{\textbackslash}Program~Files}以下に -インストールすると，アクセス制御の -関係でパッチなどが当てられなくなるので，\texttt{C:{\textbackslash}Program~Files}以外への -インストールが推奨されている -\footnote{\url{http://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-winbugs/#install}}。 -インストーラーは64ビットWindowsには非対応である。 - -\paragraph{Macへのインストール} - -Macでは，\textsf{Wine}を利用して\textsf{WinBUGS}を使うことができる。 -macOS用の\textsf{Wine}\footnote{\url{https://wiki.winehq.org/MacOSX}} -を利用するか，\textsf{Homebrew}\footnote{\url{http://brew.sh/}}, \textsf{MacPors}\footnote{\url{http://www.macports.org/}}などの -パッケージ管理システムを利用してインストールするのが簡単であろう。 -詳細は，ネット上の資料 -\footnote{\url{http://nhkuma.blogspot.jp/2012/12/macosx106-107winbugswiner2winbugs.html}\\ -\url{http://nhkuma.blogspot.jp/2012/12/macosx108-mountain-lion-winbugs-wine.html}}を参照するとよい。 - - -\textsf{Wine}をインストールしたら， -これを使用してWindows用インストーラー -(\texttt{WinBUGS14.exe})を起動 -し，\textsf{Wine}環境へ -\textsf{WinBUGS}をインストールする -（デフォルトでは\texttt{\textasciitilde/.wine/drive\_c/Program Files}以下にインストールされる）。 - - -\paragraph{Linuxへのインストール} -\textsf{Wine}を使用する。Ubuntuなどの主要なディストリビューションではパッ -ケージ化されているので，それを利用するのが簡単だろう。 -\textsf{Wine}を使用してインストーラー(\texttt{WinBUGS14.exe})を起動 -し，\textsf{WinBUGS}をインストールする。 - -BSDその他UNIXも基本的にはLinuxに準じる。 +%\subsubsection*{インストール} +% +%ウェブサイトからインストーラーをダウンロード +%できる。いずれの環境でも\textsf{WinBUGS}のインストール後には1.4.3パッチをあ +%て，``key''をインストールすること。``key''をインストールすることにより， +%全機能が使えるようになる。 +% +%\paragraph{Windowsへのインストール} +% +%インストーラー（\texttt{WinBUGS14.exe}）を起動して，あとは指示に従ってい +%けばインストールされる。 +%ただし，Windows Vistaでは，\texttt{C:{\textbackslash}Program~Files}以下に +%インストールすると，アクセス制御の +%関係でパッチなどが当てられなくなるので，\texttt{C:{\textbackslash}Program~Files}以外への +%インストールが推奨されている +%\footnote{\url{http://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-winbugs/#install}}。 +%インストーラーは64ビットWindowsには非対応である。 +% +%\paragraph{Macへのインストール} +% +%Macでは，\textsf{Wine}を利用して\textsf{WinBUGS}を使うことができる。 +%macOS用の\textsf{Wine}\footnote{\url{https://wiki.winehq.org/MacOSX}} +%を利用するか，\textsf{Homebrew}\footnote{\url{http://brew.sh/}}, \textsf{MacPors}\footnote{\url{http://www.macports.org/}}などの +%パッケージ管理システムを利用してインストールするのが簡単であろう。 +%詳細は，ネット上の資料 +%\footnote{\url{http://nhkuma.blogspot.jp/2012/12/macosx106-107winbugswiner2winbugs.html}\\ +%\url{http://nhkuma.blogspot.jp/2012/12/macosx108-mountain-lion-winbugs-wine.html}}を参照するとよい。 +% +% +%\textsf{Wine}をインストールしたら， +%これを使用してWindows用インストーラー +%(\texttt{WinBUGS14.exe})を起動 +%し，\textsf{Wine}環境へ +%\textsf{WinBUGS}をインストールする +%（デフォルトでは\texttt{\textasciitilde/.wine/drive\_c/Program Files}以下にインストールされる）。 +% +% +%\paragraph{Linuxへのインストール} +%\textsf{Wine}を使用する。Ubuntuなどの主要なディストリビューションではパッ +%ケージ化されているので，それを利用するのが簡単だろう。 +%\textsf{Wine}を使用してインストーラー(\texttt{WinBUGS14.exe})を起動 +%し，\textsf{WinBUGS}をインストールする。 +% +%BSDその他UNIXも基本的にはLinuxに準じる。 \subsection{OpenBUGS} @@ -282,7 +284,8 @@ \subsubsection*{インストール} 開発元ではWindows XPと8とで動作確認している。 \paragraph{macOSへのインストール} -\textsf{WinBUGS}と同様に，\textsf{Wine}を導入して，Windows版をインストールする。 +\textsf{Wine}を導入して，Windows版をインストールする。 +具体的な方法は，「Mac OS X 上で OpenBUGS+R2OpenBUGS を使用する - MMAye's blog」(\url{http://mmays.hatenablog.com/entry/2013/11/20/232330})などを参照。 \paragraph{Linuxへのインストール} Linux用ソースパッケージを展開して，コンパイル・インストールする。 @@ -353,7 +356,7 @@ \subsection{Stan} \subsubsection*{インストール} RStanはCRANからインストール可能。 -開発環境（C++など）が必要となる。 +モデルをコンパイルして実行するので，開発環境（C++など）が別に必要となる。 詳細はウェブサイトやマニュアルを参照。 @@ -471,16 +474,16 @@ \subsubsection{JAGSによるポアソン回帰} \end{lstlisting} % GLMでは1行で済んだものが，\textsf{BUGS言語}ではかなり長くなってしまった。 -しかし，このような記法により，複雑なモデルも記述できる。 +しかし，このような記法により，複雑なモデルも記述できるのである。 -1行目の``\texttt{model\{}''と，9行目の``\texttt{\}}''で囲まれた部分がモデルの定義となる。 -まず，3〜5行目が尤度の部分となる。 +1行目の``\texttt{model \{}''と，9行目の``\texttt{\}}''とで囲まれた部分がモデルの定義となる。 +そのうち，3〜5行目が尤度の部分である。 \texttt{for}ループの構文で，データ\texttt{X[i]}の添え字\texttt{i}の範囲が $i \in 1, 2, \dots, N$であることを示している。 \texttt{dpois}はポアソン分布であり，\texttt{dpois(lambda)}は，平均\texttt{lambda}の ポアソン分布ということになる。 ``\texttt{\textasciitilde}''(チルダ)は，左辺の変数が右辺の確率分布に従うことをしめす。 -すなわち，4行目は，\texttt{X}が，平均\texttt{lambda}のポアソン分布に従うことを表している。 +すなわち，4行目は，\texttt{X[i]}が，平均\texttt{lambda}のポアソン分布に従うことを表している。 8行目で，パラメーター\texttt{lambda}の事前分布を定義している。 \texttt{dunif}は一様分布であり，ここでは，0〜1000という範囲の一様分布を @@ -533,14 +536,17 @@ \subsubsection{JAGSによるポアソン回帰} Loaded modules: basemod,bugs \end{lstlisting} -続いて，初期値を定義する。 +続いて，MCMC計算の初期値を定義する。 これらは省略可能で，省略した場合には自動的に生成される。 -ただし，自動的に生成された値のせいでその後の計算がうまくいかないこともある。 -この例題では，初期値についてさらに，擬似乱数のタネ(\texttt{.RNG.seed})と +ただし，自動的に生成された値のせいでその後の計算がうまくいかないということもあるので， +注意が必要である。 +この例題では，パラメーター(\texttt{lambda})の初期値のほか， +擬似乱数のタネ(\texttt{.RNG.seed})と 擬似乱数発生アルゴリズム(\texttt{.RNG.name})も明示的に指定するようにしている。 -これらも省略可能である。 -しかし，パラメータの初期値と擬似乱数を固定しておくことで，再生可能(reproducible)， +これらも省略可能であるが， +パラメーターの初期値と擬似乱数を固定しておくことで，再生可能(reproducible)， すなわち同じ計算を繰り返したときに同じ結果が得られるようになる。 +% \begin{lstlisting} > inits <- vector("list", 3) > inits[[1]] <- list(lambda = 1, @@ -590,7 +596,6 @@ \subsubsection{JAGSによるポアソン回帰} % 引数の\texttt{variable.names}で，結果を保存するパラメータを指定する（ここでは\texttt{lambda}のみ）。 \texttt{n.iter}で，サンプルを取得するマルコフ連鎖の長さを指定する。ここでは50としている。 -\texttt{variable.names}引数で，サンプルを取得するパラメーターを指定する。 このほか，サンプリング間隔として\texttt{thin}引数があるが，ここではデフォルトの1を使用する（すなわちすべてのサンプルを取得する）ので，省略している。 @@ -599,7 +604,7 @@ \subsubsection{JAGSによるポアソン回帰} > traceplot(post1, col = 1, las = 1) \end{lstlisting} % -結果は図\ref{fig:trace1}。 +結果は図\ref{fig:trace1}である。 \begin{figure}[hbtp] \begin{center} @@ -609,7 +614,7 @@ \subsubsection{JAGSによるポアソン回帰} \label{fig:trace1} \end{figure} -別の初期値で試してみよう。 +別の初期値と擬似乱数でも試してみよう。 \begin{lstlisting} > inits[[2]] <- list(lambda = 30, + .RNG.seed = 2, @@ -619,7 +624,7 @@ \subsubsection{JAGSによるポアソン回帰} + .RNG.name = "base::Mersenne-Twister") \end{lstlisting} \texttt{lambda}の初期値は2番目では30，3番目では100としている。 -また，乱数のタネは2番目では2，3番目では3としている。 +また，擬似乱数のタネは2番目では2，3番目では3としている。 この設定で\textsf{JAGS}を実行する。 \texttt{jags.model()}関数の\texttt{n.chains}引数の値を3としている。 @@ -644,9 +649,9 @@ \subsubsection{JAGSによるポアソン回帰} \end{lstlisting} -軌跡を表示させる。結果は図\ref{fig:trace3}。 黒実線が1番目の初期値，赤点線が2番目の初期値，緑点線が3番目の初期値のものである。 +軌跡を表示させる。結果は図\ref{fig:trace3}。 黒実線が1番目の初期値，赤点線が2番目の初期値，青点線が3番目の初期値のものである。 \begin{lstlisting} -> traceplot(post2, las = 1) +> traceplot(post2, las = 1, col = c(1, 2, 4)) \end{lstlisting} \begin{figure}[hbtp] @@ -717,11 +722,11 @@ \subsubsection{JAGSによるポアソン回帰} 図\ref{fig:plot}をみると， 3本の連鎖がよく混ざりあっているのがわかる。このように，初期値や擬似乱数系列に -よらず，各連鎖がよく混ざりあっているのは，MCMC計算がうまくいったことを示している。 +よらず，各連鎖がよく混ざりあっているのは，MCMC計算がうまくいった（定常分布に収束したマルコフ連鎖から事後分布のサンプルを取得できた）ことを示している。 数値的に収束の診断をおこなう方法は後述する。 %ページ調整 -\pagebreak +%\pagebreak \texttt{summary()}関数で結果の要約を表示させる。 \begin{lstlisting} @@ -747,11 +752,10 @@ \subsubsection{JAGSによるポアソン回帰} % ベイズ推定ではパラメーターの推定値も確率的に（事後分布として）与えられる。 その代表値としては一般に平均や中央値が用いられる\footnote{このほか，事後分布のモードを用いる場合もある。これをMAP(Maximum A Posterior)推定値という。}。 -また，95\%信用区間（あるいは「ベイズ信頼区間」）\footnote{「頻度主義」統計における信頼区間とは実は違うもの。}も同時に使われることが多い。 +また，95\%信用区間（「ベイズ信頼区間」や「確信区間」ともいう）\footnote{「頻度主義」統計における「信頼区間」とは実は意味するものが違う。}も同時に使われることが多い。 この例では， $\lambda$の事後分布の平均値（事後平均）が3.84，中央値が3.82，95\%信用区間が -3.00〜4.81と推定された。事後平均に対応する$\lambda$の値は$\exp(1.33)=3.78$と -いうことになる。 +3.00〜4.81と推定された。 事後平均・中央値とも，GLMによる最尤推定値とだいたい一致した。 \paragraph{収束診断} @@ -782,13 +786,13 @@ \subsubsection{JAGSによるポアソン回帰} lambda 1 1 \end{lstlisting} -$\Hat{R}$値は1であった。図\ref{fig:plot}の結果とあわせ， +$\Hat{R}$の値は1であった。図\ref{fig:plot}の結果とあわせ， うまく収束したと言えそうである。 \paragraph{情報事前分布} 次に，情報を持った事前分布（情報事前分布）を使った例をみてみよう。 -今度は，$\lambda$が大きな値（10以上とか）は取らないだろうと事前の知識で +今度は，$\lambda$の値が比較的小さいと事前の知識で わかっていたとする。 そこで，$\lambda$の事前分布として$\mathrm{Gamma}(2, 2)$を指定した。 この分布の平均値は1，分散は0.5である。 @@ -866,8 +870,8 @@ \subsubsection{JAGSによるポアソン回帰} \end{lstlisting} -図\ref{prior_posterior}は，事前分布（上の点線）が尤度（下の実線）で更新されて -事後分布（上の実線）となることを示している。 +図\ref{prior_posterior}は，事前分布（上図の点線）が尤度（下図の実線）で更新されて +事後分布（上図の実線）となることを示している。 \begin{figure}[htbp] \begin{center} @@ -1131,9 +1135,6 @@ \subsubsection{Stanによる例} > library(rstan) \end{lstlisting} -%ページ調整 -\vspace{3zw} - モデルの定義。 ここでは，\texttt{example2.stan}というファイルでモデルを記述している。 % @@ -1151,9 +1152,7 @@ \subsubsection{Stanによる例} } transformed parameters { - vector[N] logit_p; - - logit_p = beta + beta_x * X; + vector[N] logit_p = beta + beta_x * X; } model { @@ -1176,9 +1175,8 @@ \subsubsection{Stanによる例} また，行末には必ずセミコロン(;)をつける。 %Stan 2.10からは，代入には``\texttt{<-}''ではなく``\texttt{=}''を使用することが推奨されるようになった。 -21行目の\texttt{binomial\_logit}分布は，logitスケール($-\infty$〜$\infty$)のパラメーターを引数にとる2項分布である。 -ここで使用するパラメーター\texttt{logit\_p}は，\texttt{transformed parameters}ブロック中の16行目で -定義されている。 +19行目の\texttt{binomial\_logit}分布は，logitスケール($-\infty$〜$\infty$)のパラメーターを引数にとる2項分布である。 +ここで使用するパラメーター\texttt{logit\_p}は，\texttt{transformed parameters}ブロック中の14行目で定義されている。 \textsf{R}コードで，連鎖の数，初期値，保存するパラメーターを設定する。 \begin{lstlisting} @@ -1190,7 +1188,7 @@ \subsubsection{Stanによる例} > pars <- c("beta", "beta_x") \end{lstlisting} -\texttt{stan()}関数によりStanを実行して，結果を\texttt{fit}というオブジェクトに格納する。 +\texttt{stan()}関数により\textsf{Stan}を実行して，結果を\texttt{fit}というオブジェクトに格納する。 \begin{lstlisting} > fit <- stan(model_code = example2_code, + data = list(X = x, Y = y, N = n, K = k), @@ -1242,7 +1240,7 @@ \subsubsection{うまくいかないとき} \item 初期値がおかしくないか? \begin{itemize} \item 例: ガンマ分布なのに負の値を与えている。 -\item 数値的にあり得ない値というわけでなくても，あまり極端な値だとエラーが発生することがある。 +\item 数学的にあり得ない値というわけでなくても，あまり極端な値だとエラーが発生することがある。 \end{itemize} \item モデルに誤りがないか? \begin{itemize} @@ -1263,7 +1261,7 @@ \subsubsection{うまくいかないとき} \item モデルが必要以上に複雑でないか? \item 余分なパラメーターがないか? \end{itemize} -\item マルコフ連鎖の長さを長くして，サンプリング間隔を長くする。 +\item マルコフ連鎖の長さを長くする。 ただし必然的に計算時間は長くなる。モデルの改良のような本質的な 改善ではないが，もう少し収束をよくしたいというようなときには 有効なこともある。 @@ -1292,7 +1290,7 @@ \subsection{線形混合モデル，あるいは階層ベイズモデル} \end{breakbox} \end{minipage} -\vspace{1zw} +\vspace{2zw} \paragraph{データ} 今回のデータはあらかじめファイルに保存してある。``\texttt{example3.csv}''から読み込む。 @@ -1300,9 +1298,6 @@ \subsection{線形混合モデル，あるいは階層ベイズモデル} > data <- read.csv("example3.csv") \end{lstlisting} -%ページ調整 -\vspace{2zw} - 最初の10行の内容を表示してみると，以下のとおり。 %\vspace{0.5em} \begin{lstlisting} @@ -1339,9 +1334,6 @@ \subsection{線形混合モデル，あるいは階層ベイズモデル} \textsf{JAGS}を使用してパラメーター推定をおこなうこととする。 使用する\textsf{R}スクリプトは``\texttt{example3.R}''である。 -%ページ調整 -%\vspace{3zw} - まず，データを行列の形に変換する（BUGSのモデルで使用する形式にあわせるため）。 各行が，各ブロックのデータとなる。 @@ -1355,7 +1347,7 @@ \subsection{線形混合モデル，あるいは階層ベイズモデル} \end{lstlisting} %ページ調整 -%\vspace{3zw} +\vspace{4zw} \texttt{x1}の内容を確認する。 @@ -1431,7 +1423,7 @@ \subsection{線形混合モデル，あるいは階層ベイズモデル} このモデルでは， ブロックによるランダム効果\texttt{e.B[]}の事前分布は，平均が0，分散が1/\texttt{tau.B} ($\texttt{tau.B}=1/\texttt{sigma.B}^{2}$)の -正規分布としている（18行目）。 +正規分布としている（19行目）。 また，\texttt{sigma.B}の事前分布は[0, 10000]の一様分布としている（28行目）。 すなわちここでは，\texttt{tau.B}(\texttt{sigma.B})が， パラメーター\texttt{e.B[]}の事前分布を決める @@ -1627,7 +1619,7 @@ \subsubsection*{Nested Indexing} モデルは以下のようになる。\texttt{X1}，\texttt{X2}，\texttt{Y}が1次元配列に なっているところ（4〜5行目），\texttt{B}という変数が導入されているところ（6行目）， -\texttt{e.B}のインデキシングが\texttt{B}を間に挟んでいるところ（16行目）などが +\texttt{e.B}のインデックスが\texttt{B}を間に挟んでいるところ（17行目）などが 前のモデルと異なっている。 \begin{lstlisting} @@ -1686,7 +1678,7 @@ \subsubsection*{中央化(centering)} \texttt{beta.1 * (X1[i] - X1.bar)}とする(\texttt{X1.bar}は\texttt{X1}の平均)。 ただし，\textsf{JAGS}で\texttt{glm}モジュールを使用した場合にはこの技法は不要である -(JAGS User's Manual\cite{JAGS} 4.6節)。 +(JAGS User's Manual 4.3 \cite{JAGS} 10.2節)。 %%