第 3 章 RStan によるべイズ統計解析 22 略称 項目 あなたは、国の政治をどれくらい信頼していますか。 1 つだけ〇を付けて 政治信頼する 問 7 ください。 現在のお宅のくらしむきを 1 年前と比べると、どうでしようか。 1 つだけ 暮らし向き改善 問 12 〇を付けてください。 この設問は、問 7 が「 1. いつも信頼している / 2. だいたい信頼している / 3. ときどきは信頼している / 4. まったく信 頼していない」問 12 が「 1. かなり良くなっている / 2. 少し良くなっている / 3. 変わらない / 4. 少し悪くなっている / 5. かなり悪くなっている」というものになっています。これを、問 7 , 12 については、回答が 1 か 2 の場合に 1 にすると いうダミー変数にします。これを読み込ませると、次のようになります。 > dataset く一 read. csv("utas—datal . ,header=T,row. name=l) > dataset 政治信頼する暮らし向き改善自民温度 0 LO 0 イ上 LD LO 0 5 っ 7 5 ( 0 2 ( 0 0 0 0 1 0 1 0 0 0 1 0 1 0 0 1 6 8 9 11 13 15 ( 略 ) 今回は、「政治信頼する」の項目と「暮らし向き改善」の項目が「自民温度」に与える影響を見てみます。モデルは次 の「 stan-test2. stan 」を使います。 data { int N int く 10wer=O, upper=l> ql CN] int く 1 ower=O , upper= 1 > q2 [N] real く 10wer=O , upper=100> y CN] parameters { real a rea1 bl rea1 b2 real く lowe て = 0 > sigma; transformed parameters { real mu[N] for ( Ⅱ in 1 : N) a + bl * ql [n] + b2 * q2 Cn] muCn] model { for (n in 1 : N) y ~ normal(mu[n] , sigma) Bayes Analysis Maniax フリ ーソフトで始めるべイズ統計解析

18 11 12 13 14 15 16 17 18 19 20 20 このデータをプロットすると右のようになります。このデータの Y を 0 .962450252 0 .958467720 0 .693674708 0 .693644477 0 .691644469 0 .592220489 0 .460733641 ー 0 .004462148 ー 0 . 106713568 ー 0 . 125764464 19 18 17 16 15 14 13 12 11 X で回帰するのが本節のモデリングです。このモデリングを、式で示し > test_lm く一 1m(YNX,data=dataset) で行ってみます。 とはいえ通常の回帰分析でもできますので、念のため通常の回帰分析 Y = 0 十 bX 十び てみます。 ( 3 ・ 1 ) 第 3 章 RStan によるべイズ統計解析 0 0 0 > test_lm Ca11 : lm(formula = Y Coefficients: (lntercept) ー 1 . 8719 15 X , data dataset) 0 . 1537 通常、回帰分析ではコマンド lm を使いますが、 RStan で回帰分析を 行う場合は、 Y が平均。十Ⅸ、分散の正規分布に従うものと仮定し 20 ( 3 ・ 2 ) ます。 data { int N ; real X [N] ; real Y[N] ; parameters{ real a; rea1 b; real く Iower=O> sigma; Bayes Analysis Maniax それに基づいて作成したモデルを作り、「 stan-testl. stan 」というファイルに保存します。内容は次の通りです。 フリーソフトで始めるべイズ統計解析

第 2 章 R によるべイズ統計解析の基本 10 > 1ibrary(LearnBayes) 次に、成功した回数と失敗した回数をベクトル形式で読み込ませ、これを事後分布を求めるコマンド pdisc に入れ ます。 > data く一 c ( 12 , 18 ) > post く一 pdisc(), prior, data) これで事後分布の予測が終わります。早速、 post の中身を見てみましよう。 > round(post,5) [ 1 ] 0 .00000 0 .17092 0 .82897 0 .00012 0 .00000 事後分布は次のような感じになりました。実験での成功率が 12 / ( 12 十 18 ) = 0.4000 だったので、 0.5 以下にさらに 偏った形になっています。 1 0.5 0.75 0.25 0 成功率が上記のものである確度 0.0000 0.1709 0.8289 0.00 ( ) 0 0.0001 事前分布と事後分布の違いを見てみるために、事後分布を出力してみます。 > plot(p,post,type="h") 2 つの分布を並べてみると、次のようになります ( 左 : 事前分布、右 : 事後分布 ) 。やはり、事後分布は 0.4 の近くに かなり偏ることが分かるはすです。 成功率 lSOd 」 0 に d 1 -0 0 8 0-6 0 4 0 、 2 0 ℃ 0 第 0 毛 0-4 0-2 2.3 べータ分布を用いる方法 先ほどの作業では事前分布が離散的な場合を使用しましたが、続いて連続的な場合を使います。この例題のような問 題を解く場合、よく使われる分布がべータ分布です。べータ分布は、次の確率密度関数で示される分布です。 0 ー 1 ( 2 ・ 1 ) ただし、 B い , の一 / ( 泳いの = ーソフトで始めるべイズ統計解析 Bayes Analysis Maniax フリ

34 第 4 章べイジアンネットワーク い / 4. どちらかと言えばよくやっているとは思わない / 5. よくやっているとは思わない / 」、問 11 が「 1. かなり良い / 2. やや良い / 3. どちらでもない / 4. やや悪い / 5. かなり悪い」、問 12 が「 1. かなり良くなっている / 2. 少し良くなってい る / 3. 変わらない / 4. 少し悪くなっている / 5. かなり悪くなっている」というものになっています。 こではこれを質的変数として使うため、 Excel を用いて、問 7 については 1 , 2 を「 1 信頼する」、 3 , 4 を「 2 信頼し ない」、問 8 は 1 , 2 を「 1 評価する」、 3 を「 2 中立」、 4 , 5 を「 3 評価しない」、問 11 , 12 は 1 , 2 を「 1 ( 景気 / 暮らし向き ) よい」、 3 を「 2 中立」、 4 , 5 を「 3 悪い」に変換します。なお、一つでも無回答があるデータは除きました。このデータ を utas-data3. csv という csv ファイルに保存します。そしてそれを datasetl という文字に読み込みます。 > datasetl く一 read. csv("utas_data3. csv" ,header=T,row. name=l) > datasetl 安倍政権 暮らし向き 1 6 8 9 11 13 政府 1 信頼する 2 信頼しない 1 信頼する 2 信頼しない 3 評価しない 3 景気悪い 3 暮らし向き悪い 2 信頼しない 3 評価しない 3 景気悪い 3 暮らし向き悪い 1 評価する 3 景気悪い 2 暮らし向き中立 1 評価する 3 景気悪い 1 暮らし向きよい 2 中立 3 景気悪い 3 暮らし向き悪い 2 信頼しない 3 評価しない 3 景気悪い 1 暮らし向きよい データの特徴を図に示します。この特徴を踏まえた上で、分析を行っていきましよう (N = 1 , 772 ) 。 合は、データセットをそのままコマンド network の中にぶち込んでしま まず、全ての組み合わせを分析対象とする場合を説明します。この場 4.3 全ての組み合わせを描画する > ggpairs (datasetl) > 1ibrary(GGa11y) > 1ibrary(ggp10t2) は前章でも用いたコマンド ggpairs で描画しました。 なお、 0 この図 えば済みます。 > netl く一 network(datasetl) > plot(netl) 0 囲みの下段にある pl 。 t コマンドで描画してみると右のようなものが 現れるはすです。この段階ではまだ矢印は描かれていません。 次に、事前分布を設定します。次のように、コマンド jointprior に 先ほど作成したネットワークを入れると、次のように動作するはすです。 > netl—prior く一 jointprior(netl) lmaginary sample size: 162 0 それでは、いよいよネットワークの描画に映ります。コマンド learn を使ってべイズ学習を行います。なお入力す る際に、最後の「 $ nw 」は描画をする際に必要になりますので、忘れないように入れてください。 Bayes AnaIysis Maniax フリ ーソフトで始めるべイズ統計解析

これを走らせた結果が下記のものです。 lnference for Stan model : stan_test5. 4 chains , each with iter=2000 ; warmup=1000 ; thin=l ; 28 42. 11 0 . 03 2 . 14 37 . 71 40 .68 42 . 16 43.56 46.30 42.80 0 . 01 0 .82 41 . 17 42 .24 42 .80 43 .36 44.35 44.89 42.91 44. 19 44. 89 45.60 46.91 22.30 0 . 07 4 . 16 14.01 19 .47 22 . 29 25.09 30.44 18 .57 15 .74 17 . 59 18 . 57 19 .54 21.40 18 .20 14.80 17 . 02 18 . 20 19 .40 21.58 12 . 53 0 . 10 6 . 13 12 . 57 16 .62 24.52 13 .46 12 . 15 13 .46 14 . 79 17 . 22 0 .03 2 . 17 6 . 53 8 . 01 9 . 47 12 . 21 0 . 13 8 . 07 ー 14 .33 ー 3 . 28 7 . 12 17 . 88 2 . 56 0 .04 2 . 66 ー 2 .61 0 .75 2 .57 4 . 32 0 .05 3 .42 7 . 99 5 .68 8 .01 10 . 33 14 .45 19 .44 0 . 01 18 .77 19 .21 19.44 19 . 67 20 . 15 + ,y=dataset$ 自民温度 ,ql=dataset$ 政治信頼する ,q2=dataset$ 景気改善 ,q3=dataset$ 暮らし向き改善 ) > dataset—stan く一 list (N=nrow(dataset) ,K=max(dataset$ 都道府県 ) ,Pref=dataset$ 都道府県 > dataset く一 read. csv("douj in066_testdata6. csvl' ,header=T,row. name=l) これを用いて次のように R に入力します。値 T を指摘するのも忘れないようにしてください。 dataset_fit く一 stan(fi1e="stan—test5. stan" , data=dataset_stan) 第 3 章 RStan によるべイズ統計解析 post—warmup draws me an S e per chain=1000 , tOtaI post—warmup draws=4000. 75 % 97 . 5 % n—eff Rhat aC1] a [ 2 ] a [ 3 ] bl [ 1 ] bl [ 2 ] bl [ 3 ] b2 [ 1 ] b2 [ 2 ] b2 [ 3 ] b3 [ 1 ] b3 [ 2 ] b3 [ 3 ] SIgma lp— me a n 0 . 03 0 . 03 0 . 02 0 . 02 sd 1 .03 1 . 45 1 . 75 1 . 94 0 .35 25 % 8 . 28 0 .45 9 .73 3 .86 1 .32 50 % 1 .82 8 . 00 1 . 89 ー 5555.02 7 . 80 0 .07 2 . 60 ー 5561.24 ー 5556.47 4000 4000 4000 3510 4000 4000 4000 4000 4000 4000 4000 4000 4000 1584 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Samp1es were drawn 11S1ng NUTS (diag-e) at Tue Ju1 25 ー 5554.70 ー 5553. 16 ー 5551 . 06 20 : 37 : 33 2017 . eff meaSUre Of effective Size, and Rhat is the potential scale reduction factor 0 Ⅱ split chains (at convergence , Rhat=l) . し、それが各地域区分のパラメータに影響を与えているという過程を置きます。このとき、モデルは次のように設定し これを基本として、階層モデルに発展させていきます。回帰分析のモデルのさらに上に、全体の平均と分散を設定 ちょっと惜しいです。 果も見てみましよう・・・・・・見た感じでは、あまり収束していないようですね・・・ 見られますが、それ以外については東部日本の他の地方とあまり変わらない値となっています。念のため MCMC の結 こうして見ると、東北地方のパラメータは、 bl ( 政治を信頼するかどうかに付随する項 ) は他の地方と若干の違いが ます。名前は「 stan-test6. stan 」とします。 int く 1 owe て = 0 int く 1 ower=0 int く 10wer=0 , int K int N data { , upper=l> q3CN] , upper=l> q2CN] upper=l> ql [N] Bayes Analysis Maniax フリ ーソフトで始めるべイズ統計解析

39 4.5 数字で示されるパラメータがある場合 物き 安権 ー■ ー 0 ■を■ー ■■ - 当 0- Co 「 0 記 1 0445 Co .0 C 「 0435 C 「 01 18- 0 18 : ] 0 18 ) : 0 ! 0 0 ・ - 刈と 0 ( ) 工 1010 0 ) ) 210 ~ ) 】 0 工 ! ・一・一価・一 : 0 : 5 印乃 180 乃門 180 : 5 乃 180 : 5 知門 18 図 4.3 サンプルデータ 2 の散布図行列 4.5 数字で示されるパラメータがある場合 こまでは項目がすべて質的なパラメータ (factor 形式 ) であるデー タを使ってきましたが、べイジアンネットワークの描画は数値で示され るパラメータ (numeric 形式 ) があってもできます。今度は、先ほどの データに加えて、問 15 で問われている、自民党、民主党、公明党、共産 党への感情温度 ( 0 ~ 100 で、 100 に近いほど強い支持、 0 に近いほど強 い不支持 ) を組み入れることにします ( ここでも 1 つでも無回答がある ものは除外。、 = 1 , 474 ) 。データは utas-data4. csv という csv ファ イルに入れます。また今回も、安倍政権を評価するかどうかを最終的な 到達点とするため、「安倍政権」から他のデータへの因果関係は計算しな いことにします。 network コマンドを使って描写を行うと、データが数 値で示されているものは黒ではなく白の丸で描かれています。

3.4 重回帰分析 政治イ言頼する 23 0 0 0 0 0 ・ 0 ・に ? 0 00 ) ・第 0 0 0 0 0 暮らし向き改善 サンプルデータの散布図行列 3 generated quantities { real y-pred CN] fO て ( Ⅱ in 1 : N) y-pred [n] = normal-rng ()u [ Ⅱ ] , 図 3.2 sigma) ; 今回は予測を行うため、モデルを一端 mu という字に格納して、その後モデルのパラメータを求めると共にモデルか dataset-stan く一 list (N=nrow(dataset) ,y=dataset$ 自民温度 ,ql=dataset$ 正攵治信頼する 民温度」を代入するために、次のように入力して、サンプリングを行います。 らの予測値 y ー pred も出すようにしています。この ql に「政治信頼する」を、 q2 に「暮らし向き改善」を、 y に「自 q2=dataset$ 暮らし向き改善 ) # $ sink("dataset—fit2. txt") dataset—fit く一 stan(fi1e="stan_test2. stanl' 結果はこのような感じです。 lnference for Stan model : stan_test2. 4 chains , each with iter=2000 ; warmup=1000 ; data=dataset_stan) thin=l ; post-warmup draws per chain=1000 , t0ta1 post-warmup draws=4000.

40 > dataset3 く一 read. csv("utas_data4. csv" ,header=T, row. name=l) > net3 く一 network(dataset3) > net3—prior く一 j0intprior(net3) lmaginary sample size : 162 第 4 章 べイジアンネットワーク > net3-ban ← as . matrix (rbind ( c ( 2 , 1 ) , c ( 2 , 3 ) , c ( 2 , 4 ) , c ( 2 , 5 ) , c ( 2 , 6 ) , c ( 2 , 7 ) , c ( 2 , 8 ) ) ) > banlist (net3) く一 net3—ban > p10t (net3) これを使って描画を行いましよう。 > net3—nw く一 1earn(net3 ,dataset3 ,net3_prior) $nw # $ > net3—post く一 autosearch(net3—nw, dataset3 ,net3_prior,trace=TRUE) Score: -30869.6 ネットワークはこのようなものになりました。 ( 略 ) Relscor 倍政 杲気 . らし向 731 e -13 、日月 : 品 主温 こうして見ると、各政党への感情温度は安倍政権への評価には直接的にも間接的にも影響しておらす、むしろ政府へ の評価や景気・暮らし向きの認識が、安倍政権への評価と各政党への感情温度の両方に影響を及ばしているという経路 4.6 付録 : 各モデルの中身 が見て取れます。 各モデルの中身は、 localprob コマンドで呼び出すことができます。下記に本書で分析したモデルの全容を示しま すので、ご参考になさってください。 Bayes Analysis Maniax フリ ーソフトで始めるべイズ統計解析

model { for (n in 1 : N) ql bernoulli(q[n] ) これを走らせてみます。 dataset—fit く一 stan(fi1e="stan—test3. stan" 結果は次の通りになりました。 lnference for Stan model : stan_test3. 3.5 ー 3 .76 0 . 11 ー 4 . 66 0 . 05 0 . 01 0 . 01 ー 4 . 04 ー 4 . 19 0 . 46 0 . 06 0 . 25 0 . 12 0 .54 0 . 36 0 . 46 0 . 10 0 . 09 0 . 01 0 . 62 0 . 01 ー 4 .35 0 . 06 0 . 11 0 .45 0 . 09 0 . 01 0 . 60 0 . 01 ー 4 . 20 0 . 01 0 . 23 0 . 01 0 . 23 0 . 00 0 . 00 0 . 00 0 . 01 0 . 00 0 . 01 0 . 00 0 . 02 0 . 00 0 . 05 0 . 00 0 . 02 0 . 00 0 . 03 0 . 00 0 . 01 0 . 00 0 . 00 0 . 00 0 . 02 0 . 00 0 . 00 0 . 03 inv-logit(a + bl * q2 [ Ⅱ ] + b2 * y ) q [ 司 for ( Ⅱ in 1 : N) real q[N] transformed parameters { real b2 real bl real a paramet ers { ロジスティック回帰分析 25 data=dataset—stan) 4 chains , each with iter=2000 ; warmup=1000 ; thin=l ; post—warmup draws per chain=1000 , total post-warmup draws=4000. 97.5 % n-eff Rhat bl b2 qC1] q [ 2 ] q[3] q[4] q [ 5 ] q[6J q[7] q[9] q [ 10 ] ( 略 ) lp- mean se_mean 0 . 02 0 . 62 0 . 02 0 . 09 0 . 10 0 .46 0 . 36 0 . 54 0 . 12 0 . 25 0 .06 0 .46 ー 823 . 83 sd 2 . 5 % 0 . 00 0 . 23 0 . 10 0 . 50 0 . 27 0 .43 0 . 06 0 . 07 0 . 58 25 % 0 . 08 0 . 33 0 . 53 0 . 24 0 . 30 ー 824.41 at Tue 50% 75 % 0 .02 0 . 63 0 . 02 0 . 10 0 . 12 0 .47 0 .40 0 . 55 0 . 13 0 . 26 0 . 06 0 . 62 0 . 91 0 .07 0 .28 0 . 14 0 .58 0 .47 0 .50 0 . 16 0 . 02 0 . 66 0 .02 1 . 23 ー 827 . 16 1259 1592 1284 2520 1430 2955 1636 4000 1520 1368 1301 2420 1301 1312 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Samp1es were drawn using NUTS (diag—e) parameter , n_eff ー 823 .53 ー 822.93 ー 822 .45 of effective Ju1 25 20 : 35 : 28 2017. and Rhat is the potential scale reduction factor 0 Ⅱ split sample Size, chains (at

4 chains , each with iter=2000 ; warmup=1000 ; thin=l ; 30 これを使って分析してみましよう。 > dataset_fit く一 stan(fiIe="stan_test6. stan" この分析の結果はこのようになりました。 lnference fO て Stan model : stan_test6. 97 . 5 % n-eff 35 . 52 10 . 90 51 . 22 42 . 55 43.79 45 . 12 57.21 43 . 20 0 . 24 39.79 42 . 18 43.26 44 . 33 45.84 42 . 96 0 . 03 41.47 43.04 42. 50 43 . 39 44.34 44.47 0 . 04 42.77 44.35 43.93 45. 03 46.43 22.35 2 . 98 21 . 19 17 .66 19 . 10 20.69 54.25 0 .24 20. 16 15 .39 17 .98 19 . 68 21.79 27.21 18 .62 0 .08 16 . 12 17 .72 18 . 56 19 . 61 21 . 23 18 .53 0 . 15 15 .38 17 .55 18 .46 19 . 59 21 .58 23. 14 96 . 61 33. 10 ー 1 . 86 9 . 01 11 . 74 17 . 11 476.39 12 .49 0 . 42 4 . 58 3 . 31 12 . 30 15 . 62 21.56 12 . 94 0 . 16 11 . 57 13 . 05 14 . 12 16 . 69 8 . 81 0 . 09 2 . 11 7 .55 8 .78 10 . 18 12 . 94 1 . 06 19 . 76 ー 21 . 25 4 .86 31 . 60 0 . 20 ー 9 . 26 0 . 00 3 . 18 14 . 35 1 . 01 0 . 12 ー 1 . 68 3 .53 1 . 01 3 . 17 0 . 37 4 .83 7 .34 12 . 99 26 . 58 0 . 31 2 .68 12 . 55 62 .80 0 . 25 2 . 62 5 . 81 80.33 51 .39 46 . 31 257 .24 0 . 69 3 .09 7 . 05 13 . 52 591 .30 17 .73 2 . 39 58 . 59 0 .82 3 .09 6 . 56 14 . 66 104.02 1 . 01 19 .39 0 . 04 0 . 33 18 .80 19 . 15 19 . 38 19 . 61 20.07 0 . 77 5 . 81 第 3 章 data=dataset—stan) RStan によるべイズ統計解析 post—warmup draws mean per chain=1000 , tOtaI post—warmup se_mean draws=4000. 50 % a0 a[l] a [ 2 ] a [ 3 ] bl ー 0 bl [ 1 ] bl [ 2 ] bl [ 3 ] b2 ー 0 b2 [ 1 ] b2 [ 2 ] b2 [ 3 ] b3 ー 0 b3 [ 1 ] b3 [ 2 ] b3 [ 3 ] S lgma_ S 1 gma_ S 1 gma— SIgma— S 1 lp- sd 1 .58 0 . 73 0 .92 2 .99 1 .32 1 . 53 1 . 92 5 . 52 2 .30 2 . 5 % 9 .62 9 . 27 4 . 63 1 . 14 ー 117 .47 25 % 75 % 24 . 28 170 . 66 bl b2 b3 4 . 37 2 . 93 3 . 20 7. 18 ー 5570 .84 9 .78 1 . 53 1 .72 1 .27 0 .94 ー 5583.25 -5574.73 ー 5570. 12 ー 5566.42 8 . 33 5 . 86 4 . 70 9 . 59 9 .35 6 . 96 20 : 40 : 07 2017. 7 . 57 246.51 ー 5560. 80 22 45 675 637 51 160 278 107 17 121 145 575 346 775 390 74 49 117 31 602 67 58 Rhat 1 .04 1 .06 1 . 10 1 .03 1 .06 1 .05 1 . 01 1 . 02 1 . 03 1 . 03 1 . 20 1 . 03 1 .02 1 . 03 1 .09 1 . 01 1 . 01 1 .07 1 . 14 Samp1es were drawn using NUTS (diag-e) at Tue Ju1 25 Ⅱ eff measure Of effective sa.mple Size , and Rhat is the potential scale reduction factor on split chains (at convergence , Rhat=l) . このように、地方による差は b2 と b3 で、東部日本と西部日本の差が出た程度ですが、東北地方は東日本の他の都道 県とあまり変わらないものに収東しました。 Bayes Analysis Maniax フリ ーソフトで始めるべイズ統計解析