Ca11 : lm(formula = Y Min Residua1s : 20 ー 1 . 61 0 . 00 0 . 01 0 . 13 0 . 18 0 . 21 0 . 42 0 . 05 1 . 32 10.92 15.97 0 .417 seconds (Tota1) 0.185 seconds (Samp1ing) E1apsed Time : 0 . 232 seconds (Warm—up) lteration: 2000 / 2000 [ 100 % ] (Samp1ing) 第 3 章 RStan によるべイズ統計解析 と、このようにべイズ推計のための MCMC ( マルコフ過程モンテカルロ法 ) が行われました ( 実際にはもっと多く の分析が行われている ) 。これによる結果を、通常の回帰分析の結果と比較してみましよう。 > summary(test—lm) data = dataset) IQ Median 0 .04038 3Q 0 . 17107 Max 0 .26980 ー 0 . 91160 ー 0 .08164 Coefficients : Est imat e (lntercept) ー 1 .87193 0 . 15368 Signif . codes: Std . 0 . 001 14 .83 0.01036 ー 15 .08 0 . 12411 Error t value Pr ( > は l) 1 . 18e ー 11 1 .55e ー 11 0 . 05 freedom 1 Residua1 standard error: 0.2672 on 18 degrees Of Mu1tip1e R—squared: 0 .9244 , Adjusted R—squared: F—statistic: 220 0 Ⅱ 1 and 18 DF, p—value : 1 .553e > dataset_fit lnference for Stan model : stan_testl . 4 chains , each with iter=2000; warmup=lOOO; thin=l ; post—warmup draws per mean S e _mean Chain sd 2 . 5 % 0 .9202 ー 11 dra s = 4000. 97 . 5 % n—eff Rhat = 1000 , total post—warmup 25 % 50 % 75 % b ー 1 .87 0 . 15 sigma 0 . 29 14 .46 lp— 0 . 00 0 .04 0 . 00 0 . 13 ー 2 . 15 ー 1 .96 ー 1 .87 ー 1 .79 0 . 15 0 . 25 13.86 0 . 15 0 . 28 14.82 0 . 16 0 . 32 15.41 1664 1600 1537 1209 1 1 1 1 Samp1es were drawn using NUTS(diag-e) at Tue Jul 25 20 : 57 : 42 2017. FO て each parameter, 取 eff crude measure Of effective sample Size , and Rhat is the potential scale reduction factor 0 Ⅱ split chains (at convergence , Rhat=1) . また、 Stan では、 MCMC がちゃんと収束したかを調べることもできます。パッケージ g cmc を使って、次のよう 率は 2.5 % 、という具合に ) 。 モデリングでは、パラメータは平均値及び実際のパラメータがある確率が示されます ( 例えば値が 2.15 以下になる確 通常の回帰分析では、各パラメータの値と標準偏差くらいしか分からないのに対して、 RStan を使ったべイズ推計の に行います。 Bayes Analysis Maniax フリ ーソフトで始めるべイズ統計解析

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.

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 フリ ーソフトで始めるべイズ統計解析

第 3 章 RStan によるべイズ統計解析 75 % 97.5 % n—eff Rhat 44.48 45.28 3670 1 22 . 10 23 . 50 3772 1 12 .48 1 9 . 98 4000 1 20 . 07 20 . 52 4000 44.48 45 . 28 3670 1 44.48 45 . 28 3670 1 66 . 07 67. 15 4000 1 54.06 56 .48 4000 1 1 66 .07 67. 15 4000 4000 1 54.06 56.48 45.28 3670 1 44.48 44.48 45.28 3670 1 1 44 .48 45 . 28 3670 45 . 28 3670 1 44 .48 24 25 % 43 . 63 20 . 71 7 .47 19 . 60 43. 63 43. 63 64.89 51 . 51 64 .89 51 . 51 43 . 63 43 . 63 43 . 63 43 . 63 2 . 5 % 42 .84 19 .38 5 . 10 19 . 16 42 .84 42 . 84 63 . 70 49 . 11 63. 70 49 . 11 42 .84 42.84 42.84 42 .84 50% 44.04 21 .39 8 .74 19 .84 44 .04 44 . 04 65 .44 52 . 80 65 .44 52 . 80 44.04 44.04 44.04 44 .04 sd mean mean 44. 06 0 . 01 0 . 62 21 . 41 0 . 02 1 .06 8 . 75 0 . 03 1 .88 19 . 84 0 . 01 0 .35 44.06 0 . 01 0 . 62 44.06 0 . 01 0 . 62 0 . 88 65 .46 0 . 01 1 . 90 52 .80 0 . 03 0 . 88 65 .46 0 . 01 52 . 80 1 .90 0 .03 44.06 0 . 01 0 .62 44.06 0 . 01 0 .62 44.06 0 . 01 0 . 62 0 . 62 44.06 0 . 01 b 1 b2 sxgma mu[l] mu [ 2 ] mu [ 3 ] mu [ 4 ] mu[5] mu [ 6 ] mu [ 7 ] mu [ 8 ] mu [ 9 ] mu[10] ロ y-pred [ 1 ] y-pred [ 2 ] y-pred [ 3 ] y-pred [ 4 ] y-pred [ 5 ] y-pred [ 6 ] y-pred [ 7 ] y-pred [ 8 ] イ亠イ亠イ上 1 一《 1 亠 -1 一 1 亠ィー一 81 . 68 4000 81.87 4000 104.56 4000 91. 19 3820 105.25 4000 90 . 67 3956 83 . 06 3971 83 .44 3872 57. 17 57 .24 78.46 66.59 79.02 66 .39 57. 16 57 .72 44. 24 43 . 72 65 . 02 53 . 04 65 .62 53.03 44.53 43 .88 29 . 96 30.57 51.73 39 .43 52 .43 39 . 61 30 . 71 30 . 27 3 . 87 5 . 34 27 . 01 12 .92 27.91 14 .20 4 .56 3 .76 0 . 31 19 . 82 19 .74 0 . 31 0 . 32 19 . 94 0 . 32 19 .93 0 . 31 19 .73 19 . 68 0 . 31 0 . 31 19 .73 0 . 32 20 . 13 43 .77 43 .76 65 . 30 52 . 83 65 . 88 53. 12 44. 16 43.82 ( 略 ) lp- 1 1974 ー 5608.20 ー 5607. 16 ー 5606.47 ー 5605 .75 0 .03 1 .39 ー 5610 . 96 ー 5607.49 NUTS (diag—e) at Tue Ju1 25 20 : 33 : 10 2017. Samp1es were drawn 1-lSIng eff iS a crude measure Of effective For each parameter , and Rhat is the potential scale reduction factor 0 Ⅱ split chains (at convergence , Rhat=l) . このように、 y ー p て ed のところにモデルから予測された値の平均値とべイズ信頼区間が表示されます。 3.5 ロジスティック回帰分析 また、 RStan では、回帰分析の目的変数が 0 から 1 の間の値の場合に用いられるロジスティック回帰分析も可能で す。今回は政府を信頼するかどうかの値を目的変数としてみましよう。モデルを次のように書き換えます。値をコマン ド inv-logit を使って文字 q に格納して、その後 ql がベルヌーイ分布に従うというモデルに入れるというものです。 ファイル名は「 stan-test3. stan 」とします。 data { int N int く 10wer=O, upper=l> ql CN] int く 10wer=O, upper=l> q2CN] rea1 く 10we て = 0 , upper=100> y Bayes Analysis 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 フリ ーソフトで始めるべイズ統計解析