杉野 勇 (SUGINO Isamu)　お茶の水女子大学・人間発達科学専攻・応用社会学コース担当

　頻度主義においては，「確率」とは，認識する主体とは無関係に，世界に事実として存在する性質を意味している。これを「客観確率」と呼んでおく。サイコロを振った時に1の目が出る確率が1/6であるというのは，誰がそれについて考えるかとは関係のない数字であると想定されている。実際に，或具体的なサイコロで1の目が出る確率が本当に1/6であるかどうかは，原理的には確認するのは難しい。サイコロが綿密に製造されていれば，たくさんの回数そのサイコロを振っていればおおよそ1/6で安定していく事が期待されるが，無限にそのサイコロを振り続けて1の目が出る確率を記録して行っても，厳密に1/6であり続ける保証はなく，ほんの僅か歪んでいる可能性が否定出来ないかも知れない。そもそもそのような無限の試行は実際には困難であるし，振っている間に摩耗して正確なサイコロで無くなっていくかもしれない。
　これに対して「主観確率」は，主体・行為者が思念の中で付与する確率である。「私はこのサイコロで1の目が出る確率は1/6だと思う／信じる」と云うのは主観確率である。「この授業で単位が取れる確率は五分五分だな」と云うのも主観確率である。世界が事実として有する性質であると想定されている客観確率は，（その確率を考える事象が時間的・空間的・物理的に同一であれば）変化しない筈であるが，同一の事象についての主観確率は簡単に変化するし，人によって値が全く異なり得る。その事に何も不合理性はない。「確率評価が変化する」と云う言い方をしても良いだろう。
　ベイズ統計学は上の主観確率を対象としていると考えられる事もあるが，必ずしも主観的であると云う部分，或は個人差を強調する必要はない。重要なのは，主体・行為者によって異なるかどうかと云う事ではなく，確率評価が（理由を持って）変化すると云う事である。例えば，以前は或要因間関連が存在する確率を10%と見積もっていたが，何らかのデータ・情報・証拠が手に入った事によりその認識が改められ，その関連の存在に対する確率評価が40%と変化する，と云うのがベイズ原理の基本的な考え方である。「事前確率が，データによって更新されて，事後確率に変化する」と表現出来るだろう。

　データ・情報によって確率が更新される例を幾つか考えてみよう。
よくある例は，病気の検査の例である。或病気に罹患する確率が10%(.10)だとしよう。その病気の検査は，90%の確率で正しく識別するとする。言い換えれば，罹患していない人に対して陽性反応を出す偽陽性の確率が10%，罹患している人に対して陰性反応を出す偽陰性の確率も10%である（偽陽性と偽陰性の確率は必ずしも等しいとは限らない）。この時，検査で陽性となった人が，その病気に罹患している確率はどれくらいだろうか？
　一般にはこの問いに対して，検査の精度である90%に近い確率を答える人が多い。しかしベイズの考え方では，（何の情報も無い状態での）事前確率は，一般的な罹患率を考える。つまり，その人が罹患している事前確率は10%である。一方，検査で陽性が出る人は，罹患していて検査で正しく検出される人と，罹患していないのに偽陽性が出る人の二種類がいる。検査で陽性が出たと云う事は，少なくともこの二つのグループのうちのどちらかに属する事は確かなので，「検査で陽性が出た時に，実際に病気に罹患している確率」は，この二つのグループのうち前者に含まれる確率，となる。
　1000人の人がいた時に，罹患者は100人，そうでない人は900人と考えられる。これが全員この検査を受けたとすると，正しく陽性反応が出るのが100人×.90＝90人，偽陽性が900人×.10＝90人，合計で180人である。この180人に含まれる人が，実際に罹患している確率は，90/180＝50%と云う事になる。事前確率が小さい場合（事前分布で比率がかなり少ない場合），事後確率は事前確率の10%よりは大幅に高くなるが，それでもまだ検査の精度90%に比べるとかなり低い値にとどまる。
　事前確率と検査の精度を与えると，更新された事後確率を出力するRスクリプトが以下である。色々な数値を設定して結果を確認してみよう。下の実行結果では上と少し異なる数値で例示している。


　科学的知識もこうした確率の更新（信念の更改）と同じ形式で考える事が出来る。
或因果説明が正しいと云う事前確率が，その後得られた証拠に照らして更新され，事後確率が高まる（低まる），と考える。
　例えば，最初，因果関係についての理論仮説Aが正しいかどうかの確率は半々であるとする（事前確率.50）。
次に，その仮説Aが正しい時に，或事象bが観察される確率が，.75であるとする。
その仮説Aが正しくない時に，同じその事象bが観察される確率が，.25であるとする。
　もし，定性的調査の結果，その事象bが観察されたとすると，仮説Aが正しい確率はどうなるであろうか。

　理論（または理論的説明）の尤もらしさが，（その理論の観察可能な含意についての）観察結果により，事前確率から事後確率に変化するのである。

　以下の例のそれぞれがどの様な研究の成果を意味しているのか考えてみよう。

　ここで詳しく説明するのはとても無理なので，詳細は専門の書籍やウェブサイトを参照して貰う事として，最低限の知識だけを例示していくにとどめる。
MCMCpackのインストール
　もし，library(MCMCpack)のところで他に足りないパッケイジがあると云った警告が出てうまく読み込まれなければ，install.packages()関数にそのパッケイジ名を指定して追加インストールして下さい。

ステューデントのt検定と二値変数による単回帰分析の確認
　次に，第6章1-2の2群の母平均の差のステューデントt検定を改めて実行しておく。

ステューデントのt検定と二値変数による単回帰分析は完全に一致する事をもう一度確認しておこう。

MCMC; Markov Chain Monte Carlo法
　ここからがベイズ統計学の概略。大胆に簡略化して説明すると，
�@ 求めたい母数の値について，出来るだけ条件を置かない事前確率分布を指定する（一様分布や，極めて分散の大きな正規分布など）。つまり，とても広範囲の様々な値を，ほぼ同じ様な（小さな）確率で取り得ると仮定する。
�A 調査で得られたデータから，様々な母数の尤度を求める。
�B ベイズの原理に従って，事前確率×尤度＝事後確率が最大になる母数を求める。
実際には，この事後確率の分布，事後確率分布を求めるのは極めて大変なので，乱数を活用したシミュレイションによって事後分布を求める方法が，コンピュータの性能の劇的な向上と共に一般化した。その代表的なものが，マルコフ連鎖モンテカルロ(MCMC)法と呼ばれる方法である。
　MCMCでは，乱数から出発してランダムに変動する一繋がりの（しばしば膨大な数の）推定値群を得る。初期値に用いた乱数の影響が及ぶので最初の方の推定値群は母数の推定値としては好ましくないが，十分な回数の変化を経た後は，母数の推定値として使用出来ると考えられる。そうした推定値が膨大な数だけ集まれば，それを母数の事後確率分布と見做す事が出来る。
　初期値の影響が残っているかどうかを確認する為に，初期値を変えて乱数の系列生成を行い，最終的に似た様な事後分布になるかどうか（収束するかどうか）を確認する。
　ここまでを，上の幸福度と性別の例で試してみよう。使用する関数は，MCMCregress( )である。burnin=は，最初の何個の数字を捨てるか，mcmc=はその系列で推定値を何個発生させるか，thin=は何個おきに推定値を残すか，seed=は乱数の種の指定である。

これは推定値を10個(繰り返しの数 mcmc=10)しか生成しておらず，連鎖（chain）は出たらめに上下し，ヒストグラムの形も滑らかではない。通常は，繰り返しの数はもっと遥かに大きな値を指定する。以下では，mcmc=100, mcmc=1000, mcmc=10000として，グラフがどの様になっていくかを示してみよう（スクリプトは省略，グラフも真ん中の factor(sex)2 のもののみ）。



　以下では，burnin=1000として最初の1000個の数字は捨ててmcmc=100000としてそれぞれ100000個の数字からなる異なる系列を3つ生成した上で10個おきの数字だけを採用し(thin=10，3系列を1つに纏めた上で図示し，且つ summary( )で計算結果の数字を確認しよう。また，収束しているか否か（きちんとした信頼できる結果が得られているか否か）を診断する関数も付記だけしておく（収束診断の結果は省略）。


切片(Intercept)，男性に対する女性の平均差(factor(sex)2)，誤差分散(sigma2)のそれぞれについて，10000個の推定値が3系列あるので全部で30000個の推定値が得られ，その事後分布の平均(Mean)や標準偏差(SD)，2.5%点や97.5%点が表示されている。
2.5%点から97.5%点の範囲は，頻度主義ではconfidence interval（信頼区間）と呼ぶが，ベイズ統計学においては，credible interval（信用区間或は確信区間）と呼ばれている。ベイズのcredible intervalの解釈は非常に素朴で分かり易く，「母数がその区間に存在する確率が95%である。」と云ったものになる。
　この結果を，上のステューデントのt検定または単回帰分析の結果と比較してみると，新しい統計学を導入した割には，結果は殆ど一致しているので拍子抜けするかも知れない。実は，初歩的な重回帰分析くらいでは頻度主義の95%信頼区間とベイズの95%確信区間はほぼ変わらない事が多い。ベイズ統計学を本当に活用する為には，もっと複雑なモデル，前提を置かないモデルなどで柔軟な分析をするのが良いかも知れない。

　基本的に目新しい内容は出て来ないが，第9章1-1の重回帰分析をベイズで分析してみよう。
　まずは，第9章の通常の(OLS)重回帰分析を確認しておく。

　この結果には後程戻って来て貰う事として，これと同じモデルを MCMCregress( )で計算しよう。



　分析結果は，OLS重回帰分析とほぼ一致する。頻度論的推定は，事前分布に前提を置かないベイズ推定とほぼ一致する。
なお，収束診断の結果は省略したが，この場合は必ずしも良好とは言えない様である。モデルを改善する必要がありそうだ。

　第12章1-3の2項ロジスティック回帰（2項ロジット）も，MCMClogit( )でほぼ同様に分析する事が出来，事前確率分布に特に条件を設定しない限り，頻度論的なロジスティック回帰分析の結果とほぼ同様となる。
　データの読み込みや変数作成の部分は省略する。#の後のコメント分は邪魔ならなくて良い。
　まずはコンソール画面の計算結果。第12章1-3の glm( , family = binomial(link = "logit")) の結果とほぼ一致する（よっていまいち有難みは感じられないが…）。ついでに何種類かある収束判定方法の一つGelman法の計算結果も。

　次は，或時点の値とそれとlagだけ離れた時点の値の相関係数（自己相関）を，1番目のチェイン(BL1)についてだけ，切片を除いて描画したものである。0に近づいていれば問題無い。

　次のグラフは，トレースプロットと密度プロットと呼ばれるもので，前者が一つの帯の様になっていれば問題無く，密度グラフがベイズ信用区間（確信区間）を視覚化している。スペース節約の為切片は省略した。


　最後は，Gelman法の収束診断についての図示である。これも切片は省略した。1.1以下で許容可能と言われている様だ。

　上にも述べた様に，rstanのインストールや活用は初級者にはややハードルが高い。下記のサイトや書籍を参照して欲しい(いずれも2017年5月7日閲覧)。
Kentaro Matsuura, 'RStan Getting Started (Japanese): RStanをはじめよう'，
松浦健太郎, 2016, 『StanとRでベイズ統計モデリング』, 東京: 共立出版.
　ここでは，RStudioは使用せずにプレーンなRで演習する。rstanのinstallまでは出来ているものとして，ごく単純な重回帰分析を行った例を示す。演習用データで，男性ダミーmaleと年齢ageと教育年数eduから個人年収incomeを説明（推定）するモデルを作成する。
　RとStanは本来は別々のアプリケイションであり，RのパッケイジrstanはRとStanを繋ぐものであると考えると良い。従って，やや面倒であるが，Stanのコードファイルと，それを呼び出すRのスクリプトファイルの2つを作成する必要がある。
　まずはStanのコードファイル(ファイル名を MRA-3IV.stan として保存しておく)。
次に，これを呼び出すRのスクリプトである（→Mac用Rのエディタでの表示画面）。
この出力結果は長大になるので，テキストファイルで別ウィンドウで見て欲しい。
上のRスクリプトのコンソール出力結果。

　重要な部分（ bayesMRA の行）を中心に，以下に掲示する（Mac版Rのコンソール）。

　古典的／頻度論的重回帰分析での分析結果は以下である（Mac版Rのコンソール）。


　更に，ほぼ条件を揃えた MCMCregress( )での分析結果も示しておく。





　いずれの方法でも結果はほぼ同一と言って良い程度であり，この程度の単純な分析ではベイズ統計学（MCMC法による計算）やStanの有用性は余り分からないが，Stanのコードファイルを見ると何となく分かる様に，複雑なモデルを自由に組む程，Stanの有効性が高まると考えられる。

　もう少しrstanの説明を補おう。重回帰分析を例示した後で二群の母平均の差の推測をするのはやや順序が逆にも思えるが，より基本的な分析に立ち戻って初歩から理解しよう。
　頻度論の検定では，二群の母分散が等しいと言えるか言えないかで場合分けをする習慣がかつてはあり，母分散が等しければStudentのt検定，母分散が等しいとは言えなければWelchのt検定を行う事が多かった。その後，等分散であるか否かを確認せずに（検定の多重性を回避して）最初からWelch検定を行う事が主流となった。
　Stanを用いた分析では，そもそも等分散を仮定する必要が無いので（仮定したければ勿論そう出来る），最初からWelchと同じアプローチを取る事になる。
　最初にStanのコードファイル（青の背景），そして次にそれを呼び出すRのスクリプトファイル（黄色の背景）を示す。
　次は，このStanコード（bayes_welch.stan）を呼び出すRのスクリプトファイルである。
　オブジェクト bayes_welch のMac上のRでの表示。その上には計算過程のコンソール出力の最後の辺りが表示されている。


　次が stan_trace( ) の出力。問題無く収束しているかどうかを目視する。chain毎に描き分けている。


　stan_histの出力は省略し，次は stan_dens( ) の出力。


　次は plot(density(data04$diff), main="平均値の差の事後分布") の結果。


　最後に，「# 二つの平均値の事後分布をヒストグラムで重ね描き」以下の描画結果。

　男性の幸福度平均も女性の幸福度平均もそれぞれバラつきはあるが，ほぼ確実に女性の幸福度の方が高い事が読み取れる。Rスクリプトの post_prob の部分で，女性の幸福度平均 マイナス 男性の幸福度平均 の事後分布が0より大きくなる確率を計算している。結果は省略しているが，99%以上となる。
　頻度論的なスチューデントのt検定，Welch検定の結果は以下の通りである。Student検定とWelch検定の違いも大きくはないが，上のrstanの結果（オブジェクト"bayes_welch"のコンソール出力）と見比べると，Welch検定の差の95%信頼区間とベイズ分析の95%信用区間の方が僅かに近くなっている。また，ベイズ推計では，二群の母標準偏差もそれぞれ2.03，1.77と具体的に推計されている。分散比にすれば1.3倍余りである。


　結果の貼り付けだけになるが，これらを参考にすれば，MCMCpackパッケイジの plot( ) と同じかそれ以上のグラフを rstan で作成出来るだろう。

　第13章1-2 マルチレヴェルモデルで例示したモデルを，rstanでも再現してみる。
以下が，このStanコードファイルを呼び出すRスクリプトである。データの読み込みやmaleダミー変数，地点変数の作成の部分は第13章1-1と同じなので省略する。
　まずは基本となる出力結果である。第13章1-2の結果と見比べて多少係数値は相違するが，95%信用区間の広さに比べれば大した違いではないと言えるだろう。


stan_trace(bayesHLM, pars=c("b0", "b1", "b2", "b3", "b13")) の描画結果。これまでの例と比べると，やや収束が悪い（どのchainも同じ様な帯になってはいない）様に見える。


stan_hist(bayesHLM, pars=c("b0", "b1", "b2", "b3", "b13")) の描画結果。


stan_dens(bayesHLM, pars=c("b0", "b1", "b2", "b3", "b13"), separate_chains=T) の描画結果。chainごとの違いがやや目につく。


stan_ac( ) の描画結果は省略するが，自己相関の低下も遅い／十分でない可能性が伺われる。