20190725nigata-test-simple
set.seed(1234)nigata data
read table & rounding to integers
dat <- read.csv("nigata.csv")
colnames(dat) <- c("city", "sosu", "tokuhyo")
data = list(N=nrow(dat), x=round(dat$tokuhyo), n=round(dat$sosu), city=as.character(dat$city))plot
plot(data$n, data$x, type="n", xlab="投票総数", ylab="得票数", family=family); text(data$n, data$x, cex=0.5)
Computing posterior distributions (Bayesian inference)
stan
library(rstan)## Loading required package: StanHeaders## Loading required package: ggplot2## rstan (Version 2.19.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores()) fit
betabinom2 = stan_model("betabinom2.stan")
fit1 <- sampling(betabinom2, data=data, chains=8, warmup=2000, iter=20000)stan diagnostics
fit1## Inference for Stan model: betabinom2.
## 8 chains, each with iter=20000; warmup=2000; thin=1;
## post-warmup draws per chain=18000, total post-warmup draws=144000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## p[1] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## p[3] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[4] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[5] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[6] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[7] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[8] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[9] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[10] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[11] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[12] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[13] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[14] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[15] 0.00 0.00 0.00 0.00 0.00 0.00 0.01
## p[16] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[17] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[18] 0.01 0.00 0.00 0.00 0.00 0.01 0.01
## p[19] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[20] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[21] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[22] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[23] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[24] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[25] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[26] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[27] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[28] 0.01 0.00 0.00 0.00 0.00 0.01 0.01
## p[29] 0.01 0.00 0.00 0.00 0.01 0.01 0.01
## p[30] 0.01 0.00 0.00 0.00 0.01 0.01 0.01
## p[31] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[32] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[33] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[34] 0.00 0.00 0.00 0.00 0.00 0.00 0.01
## p[35] 0.01 0.00 0.00 0.00 0.00 0.01 0.01
## p[36] 0.01 0.00 0.00 0.00 0.01 0.01 0.01
## p[37] 0.01 0.00 0.00 0.00 0.01 0.01 0.01
## phi 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## lambda 1123.27 0.67 280.37 646.60 923.47 1098.67 1295.38
## alpha 8.55 0.00 2.10 4.99 7.06 8.35 9.83
## beta 1114.73 0.66 278.33 641.62 916.51 1090.23 1285.56
## x_rep[1] 902.10 0.85 322.81 382.00 672.00 865.00 1090.00
## x_rep[2] 660.90 0.62 236.83 279.00 492.00 633.00 798.00
## x_rep[3] 576.47 0.54 206.85 244.00 429.00 553.00 696.00
## x_rep[4] 507.04 0.48 182.03 214.00 377.00 486.00 613.00
## x_rep[5] 414.19 0.39 149.03 174.00 308.00 397.00 501.00
## x_rep[6] 336.23 0.32 121.20 141.00 250.00 322.00 407.00
## x_rep[7] 317.34 0.30 114.40 133.00 235.00 304.00 384.00
## x_rep[8] 302.31 0.29 109.14 127.00 225.00 290.00 366.00
## x_rep[9] 275.69 0.26 99.61 115.00 205.00 264.00 334.00
## x_rep[10] 249.38 0.24 90.24 104.00 185.00 239.00 302.00
## x_rep[11] 226.05 0.21 81.96 94.00 168.00 217.00 274.00
## x_rep[12] 215.77 0.21 78.21 90.00 160.00 207.00 261.00
## x_rep[13] 211.09 0.20 76.64 88.00 156.00 202.00 256.00
## x_rep[14] 210.77 0.20 76.40 88.00 156.00 202.00 255.00
## x_rep[15] 206.62 0.20 75.05 86.00 153.00 198.00 250.00
## x_rep[16] 206.23 0.20 74.84 86.00 153.00 198.00 250.00
## x_rep[17] 182.47 0.17 66.36 76.00 135.00 175.00 221.00
## x_rep[18] 174.39 0.17 63.50 72.00 129.00 167.00 211.00
## x_rep[19] 165.75 0.16 60.46 69.00 123.00 159.00 201.00
## x_rep[20] 149.48 0.14 54.63 61.00 110.00 143.00 181.00
## x_rep[21] 146.22 0.14 53.45 60.00 108.00 140.00 177.00
## x_rep[22] 140.11 0.14 51.32 58.00 103.00 134.00 170.00
## x_rep[23] 139.25 0.13 50.96 57.00 103.00 133.00 169.00
## x_rep[24] 133.19 0.13 48.81 55.00 98.00 128.00 162.00
## x_rep[25] 127.74 0.12 46.91 52.00 94.00 122.00 155.00
## x_rep[26] 107.59 0.10 39.76 43.00 79.00 103.00 131.00
## x_rep[27] 96.19 0.09 35.63 39.00 71.00 92.00 117.00
## x_rep[28] 48.03 0.05 18.44 18.00 35.00 46.00 59.00
## x_rep[29] 42.41 0.04 16.44 16.00 31.00 41.00 52.00
## x_rep[30] 42.38 0.04 16.39 16.00 31.00 41.00 52.00
## x_rep[31] 41.38 0.04 16.09 16.00 30.00 39.00 51.00
## x_rep[32] 29.68 0.03 11.91 11.00 21.00 28.00 37.00
## x_rep[33] 28.65 0.03 11.50 10.00 20.00 27.00 35.00
## x_rep[34] 25.04 0.03 10.23 9.00 18.00 24.00 31.00
## x_rep[35] 19.14 0.02 8.10 6.00 13.00 18.00 24.00
## x_rep[36] 17.78 0.02 7.59 6.00 12.00 17.00 22.00
## x_rep[37] 1.73 0.00 1.45 0.00 1.00 1.00 3.00
## lp__ -46080.65 0.02 4.48 -46090.33 -46083.46 -46080.33 -46077.47
## 97.5% n_eff Rhat
## p[1] 0.01 251242 1
## p[2] 0.00 243459 1
## p[3] 0.01 263989 1
## p[4] 0.01 256447 1
## p[5] 0.01 258196 1
## p[6] 0.01 261860 1
## p[7] 0.01 261094 1
## p[8] 0.01 260181 1
## p[9] 0.01 264776 1
## p[10] 0.01 249604 1
## p[11] 0.01 257052 1
## p[12] 0.01 253736 1
## p[13] 0.01 266589 1
## p[14] 0.01 257041 1
## p[15] 0.01 259572 1
## p[16] 0.01 254716 1
## p[17] 0.01 251930 1
## p[18] 0.01 258953 1
## p[19] 0.01 256887 1
## p[20] 0.01 256434 1
## p[21] 0.01 252297 1
## p[22] 0.01 258417 1
## p[23] 0.01 242746 1
## p[24] 0.01 261879 1
## p[25] 0.01 259031 1
## p[26] 0.01 263776 1
## p[27] 0.01 252740 1
## p[28] 0.01 247501 1
## p[29] 0.01 257782 1
## p[30] 0.01 256643 1
## p[31] 0.01 258634 1
## p[32] 0.01 249935 1
## p[33] 0.01 256786 1
## p[34] 0.01 235431 1
## p[35] 0.01 235602 1
## p[36] 0.01 265218 1
## p[37] 0.01 241934 1
## phi 0.01 190093 1
## lambda 1740.10 176397 1
## alpha 13.18 181941 1
## beta 1727.18 176368 1
## x_rep[1] 1635.00 145779 1
## x_rep[2] 1200.00 145293 1
## x_rep[3] 1046.00 145797 1
## x_rep[4] 921.00 145536 1
## x_rep[5] 753.00 145095 1
## x_rep[6] 612.00 145330 1
## x_rep[7] 577.00 144995 1
## x_rep[8] 551.00 145779 1
## x_rep[9] 502.00 145650 1
## x_rep[10] 454.00 145168 1
## x_rep[11] 412.00 145743 1
## x_rep[12] 394.00 145494 1
## x_rep[13] 386.00 145399 1
## x_rep[14] 384.00 145247 1
## x_rep[15] 377.00 145173 1
## x_rep[16] 377.00 145553 1
## x_rep[17] 333.00 145725 1
## x_rep[18] 319.00 145149 1
## x_rep[19] 303.00 145542 1
## x_rep[20] 274.00 145677 1
## x_rep[21] 267.00 144457 1
## x_rep[22] 257.00 144320 1
## x_rep[23] 255.00 145595 1
## x_rep[24] 244.00 145015 1
## x_rep[25] 234.00 145353 1
## x_rep[26] 198.00 145091 1
## x_rep[27] 177.00 145488 1
## x_rep[28] 90.00 144234 1
## x_rep[29] 80.00 144661 1
## x_rep[30] 80.00 145197 1
## x_rep[31] 78.00 145724 1
## x_rep[32] 57.00 144688 1
## x_rep[33] 55.00 143997 1
## x_rep[34] 48.00 143840 1
## x_rep[35] 38.00 145349 1
## x_rep[36] 35.00 145465 1
## x_rep[37] 5.00 142802 1
## lp__ -46072.84 55503 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 25 22:08:54 2019.
## 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).stan_trace(fit1, pars=c("phi", "lambda"))
stan_dens(fit1, pars=c("phi", "lambda"), separate_chains = TRUE)
#stan_hist(fit1)stan_dens(fit1, pars=c("x_rep[1]","x_rep[2]","x_rep[3]","x_rep[4]"))
stan summary
ss = summary(fit1)$summary
ss## mean se_mean sd 2.5%
## p[1] 9.027564e-03 5.468381e-07 2.740972e-04 8.498803e-03
## p[2] 1.033630e-03 2.238705e-07 1.104612e-04 8.286858e-04
## p[3] 1.095313e-02 7.333778e-07 3.768087e-04 1.022731e-02
## p[4] 1.090355e-02 7.907186e-07 4.004245e-04 1.013024e-02
## p[5] 9.552117e-03 8.171308e-07 4.152083e-04 8.760113e-03
## p[6] 9.631574e-03 8.996663e-07 4.603794e-04 8.749744e-03
## p[7] 8.477437e-03 8.687129e-07 4.438894e-04 7.626780e-03
## p[8] 7.365070e-03 8.287809e-07 4.227439e-04 6.559868e-03
## p[9] 9.076986e-03 9.564051e-07 4.921312e-04 8.141175e-03
## p[10] 9.436361e-03 1.056550e-06 5.278566e-04 8.432813e-03
## p[11] 9.885800e-03 1.110815e-06 5.631863e-04 8.812125e-03
## p[12] 6.729767e-03 9.485803e-07 4.778211e-04 5.827885e-03
## p[13] 7.880343e-03 1.014447e-06 5.237824e-04 6.888172e-03
## p[14] 6.327830e-03 9.248105e-07 4.688716e-04 5.443739e-03
## p[15] 4.885544e-03 8.187505e-07 4.171388e-04 4.101306e-03
## p[16] 8.556493e-03 1.085539e-06 5.478644e-04 7.517749e-03
## p[17] 9.143560e-03 1.196634e-06 6.006221e-04 7.997626e-03
## p[18] 5.033790e-03 8.996820e-07 4.578251e-04 4.171346e-03
## p[19] 7.916062e-03 1.154586e-06 5.851909e-04 6.808004e-03
## p[20] 7.762350e-03 1.212096e-06 6.137976e-04 6.604233e-03
## p[21] 6.641148e-03 1.133440e-06 5.693181e-04 5.571843e-03
## p[22] 9.587760e-03 1.379238e-06 7.011318e-04 8.261239e-03
## p[23] 7.469914e-03 1.255569e-06 6.186094e-04 6.310511e-03
## p[24] 1.000293e-02 1.432942e-06 7.332959e-04 8.612981e-03
## p[25] 7.094744e-03 1.237796e-06 6.299778e-04 5.915820e-03
## p[26] 7.472290e-03 1.359386e-06 6.981686e-04 6.166681e-03
## p[27] 8.063841e-03 1.522513e-06 7.654159e-04 6.633074e-03
## p[28] 5.201162e-03 1.687245e-06 8.393952e-04 3.685080e-03
## p[29] 6.366600e-03 1.920907e-06 9.752882e-04 4.596393e-03
## p[30] 5.930339e-03 1.860641e-06 9.425989e-04 4.226009e-03
## p[31] 7.582476e-03 2.115777e-06 1.076000e-03 5.615465e-03
## p[32] 8.100374e-03 2.551096e-06 1.275381e-03 5.806722e-03
## p[33] 8.536778e-03 2.617384e-06 1.326335e-03 6.139914e-03
## p[34] 4.664977e-03 2.163547e-06 1.049780e-03 2.821586e-03
## p[35] 5.375918e-03 2.542852e-06 1.234271e-03 3.201770e-03
## p[36] 6.817960e-03 2.741277e-06 1.411738e-03 4.329337e-03
## p[37] 7.070688e-03 4.819623e-06 2.370618e-03 3.177532e-03
## phi 7.631260e-03 1.063055e-06 4.634877e-04 6.762524e-03
## lambda 1.123273e+03 6.675457e-01 2.803666e+02 6.465952e+02
## alpha 8.546168e+00 4.913963e-03 2.096027e+00 4.992686e+00
## beta 1.114727e+03 6.627510e-01 2.783301e+02 6.416213e+02
## x_rep[1] 9.021016e+02 8.454746e-01 3.228108e+02 3.820000e+02
## x_rep[2] 6.608983e+02 6.213286e-01 2.368334e+02 2.790000e+02
## x_rep[3] 5.764729e+02 5.417403e-01 2.068546e+02 2.440000e+02
## x_rep[4] 5.070435e+02 4.771633e-01 1.820340e+02 2.140000e+02
## x_rep[5] 4.141929e+02 3.912348e-01 1.490265e+02 1.740000e+02
## x_rep[6] 3.362312e+02 3.179248e-01 1.211999e+02 1.410000e+02
## x_rep[7] 3.173394e+02 3.004422e-01 1.144031e+02 1.330000e+02
## x_rep[8] 3.023051e+02 2.858376e-01 1.091357e+02 1.270000e+02
## x_rep[9] 2.756924e+02 2.609951e-01 9.960632e+01 1.150000e+02
## x_rep[10] 2.493838e+02 2.368577e-01 9.024489e+01 1.040000e+02
## x_rep[11] 2.260497e+02 2.147000e-01 8.196460e+01 9.400000e+01
## x_rep[12] 2.157695e+02 2.050519e-01 7.821426e+01 9.000000e+01
## x_rep[13] 2.110930e+02 2.009914e-01 7.664050e+01 8.800000e+01
## x_rep[14] 2.107725e+02 2.004655e-01 7.640004e+01 8.800000e+01
## x_rep[15] 2.066183e+02 1.969790e-01 7.505221e+01 8.600000e+01
## x_rep[16] 2.062263e+02 1.961551e-01 7.483604e+01 8.600000e+01
## x_rep[17] 1.824659e+02 1.738460e-01 6.636375e+01 7.600000e+01
## x_rep[18] 1.743938e+02 1.666799e-01 6.350245e+01 7.200000e+01
## x_rep[19] 1.657495e+02 1.584749e-01 6.045809e+01 6.900000e+01
## x_rep[20] 1.494791e+02 1.431237e-01 5.462699e+01 6.100000e+01
## x_rep[21] 1.462186e+02 1.406398e-01 5.345370e+01 6.000000e+01
## x_rep[22] 1.401096e+02 1.350927e-01 5.132103e+01 5.800000e+01
## x_rep[23] 1.392528e+02 1.335584e-01 5.096184e+01 5.700000e+01
## x_rep[24] 1.331905e+02 1.281643e-01 4.880610e+01 5.500000e+01
## x_rep[25] 1.277430e+02 1.230301e-01 4.690543e+01 5.200000e+01
## x_rep[26] 1.075862e+02 1.043916e-01 3.976363e+01 4.300000e+01
## x_rep[27] 9.619394e+01 9.342279e-02 3.563415e+01 3.900000e+01
## x_rep[28] 4.803183e+01 4.855563e-02 1.844053e+01 1.800000e+01
## x_rep[29] 4.241133e+01 4.323340e-02 1.644352e+01 1.600000e+01
## x_rep[30] 4.237845e+01 4.302002e-02 1.639268e+01 1.600000e+01
## x_rep[31] 4.137953e+01 4.215326e-02 1.609153e+01 1.600000e+01
## x_rep[32] 2.968418e+01 3.131338e-02 1.191097e+01 1.100000e+01
## x_rep[33] 2.865074e+01 3.029935e-02 1.149767e+01 1.000000e+01
## x_rep[34] 2.503776e+01 2.697837e-02 1.023187e+01 9.000000e+00
## x_rep[35] 1.913653e+01 2.123409e-02 8.095421e+00 6.000000e+00
## x_rep[36] 1.777644e+01 1.990433e-02 7.591485e+00 6.000000e+00
## x_rep[37] 1.726493e+00 3.831794e-03 1.448004e+00 0.000000e+00
## lp__ -4.608065e+04 1.901816e-02 4.480501e+00 -4.609033e+04
## 25% 50% 75% 97.5%
## p[1] 8.840896e-03 9.025209e-03 9.211096e-03 9.570870e-03
## p[2] 9.569742e-04 1.029850e-03 1.106287e-03 1.261326e-03
## p[3] 1.069730e-02 1.094777e-02 1.120490e-02 1.170646e-02
## p[4] 1.063061e-02 1.089816e-02 1.117126e-02 1.170371e-02
## p[5] 9.268517e-03 9.545596e-03 9.828279e-03 1.038575e-02
## p[6] 9.315706e-03 9.624361e-03 9.939958e-03 1.055274e-02
## p[7] 8.173365e-03 8.470785e-03 8.773496e-03 9.368297e-03
## p[8] 7.075830e-03 7.357146e-03 7.644780e-03 8.217058e-03
## p[9] 8.739472e-03 9.068260e-03 9.404751e-03 1.006662e-02
## p[10] 9.075336e-03 9.427192e-03 9.786614e-03 1.049548e-02
## p[11] 9.501760e-03 9.874718e-03 1.025863e-02 1.102480e-02
## p[12] 6.400070e-03 6.716746e-03 7.046807e-03 7.697257e-03
## p[13] 7.519539e-03 7.868192e-03 8.226470e-03 8.942855e-03
## p[14] 6.004131e-03 6.315858e-03 6.639860e-03 7.277641e-03
## p[15] 4.597766e-03 4.874122e-03 5.161349e-03 5.734882e-03
## p[16] 8.179026e-03 8.545970e-03 8.919562e-03 9.662947e-03
## p[17] 8.732407e-03 9.132230e-03 9.538594e-03 1.035603e-02
## p[18] 4.718800e-03 5.020803e-03 5.334278e-03 5.966861e-03
## p[19] 7.514555e-03 7.901319e-03 8.302087e-03 9.104900e-03
## p[20] 7.339878e-03 7.746485e-03 8.167819e-03 9.006960e-03
## p[21] 6.248334e-03 6.626130e-03 7.016314e-03 7.798294e-03
## p[22] 9.109610e-03 9.569688e-03 1.004914e-02 1.101779e-02
## p[23] 7.043000e-03 7.450912e-03 7.877145e-03 8.729048e-03
## p[24] 9.504729e-03 9.984264e-03 1.048205e-02 1.149669e-02
## p[25] 6.659229e-03 7.073298e-03 7.512386e-03 8.380621e-03
## p[26] 6.990965e-03 7.449284e-03 7.930903e-03 8.907550e-03
## p[27] 7.532579e-03 8.040620e-03 8.567259e-03 9.630023e-03
## p[28] 4.613900e-03 5.156569e-03 5.743177e-03 6.965960e-03
## p[29] 5.685748e-03 6.319473e-03 6.992864e-03 8.419186e-03
## p[30] 5.268569e-03 5.880937e-03 6.541040e-03 7.910182e-03
## p[31] 6.832811e-03 7.531173e-03 8.277488e-03 9.835416e-03
## p[32] 7.204491e-03 8.030386e-03 8.920667e-03 1.080036e-02
## p[33] 7.607449e-03 8.471971e-03 9.392542e-03 1.130853e-02
## p[34] 3.924112e-03 4.593235e-03 5.329346e-03 6.913748e-03
## p[35] 4.503980e-03 5.292252e-03 6.150824e-03 8.035983e-03
## p[36] 5.821566e-03 6.722974e-03 7.708070e-03 9.835863e-03
## p[37] 5.375513e-03 6.823650e-03 8.478326e-03 1.242472e-02
## phi 7.318054e-03 7.615209e-03 7.926142e-03 8.594684e-03
## lambda 9.234686e+02 1.098674e+03 1.295378e+03 1.740101e+03
## alpha 7.055235e+00 8.351488e+00 9.827840e+00 1.318278e+01
## beta 9.165092e+02 1.090226e+03 1.285563e+03 1.727184e+03
## x_rep[1] 6.720000e+02 8.650000e+02 1.090000e+03 1.635000e+03
## x_rep[2] 4.920000e+02 6.330000e+02 7.980000e+02 1.200000e+03
## x_rep[3] 4.290000e+02 5.530000e+02 6.960000e+02 1.046000e+03
## x_rep[4] 3.770000e+02 4.860000e+02 6.130000e+02 9.210000e+02
## x_rep[5] 3.080000e+02 3.970000e+02 5.010000e+02 7.530000e+02
## x_rep[6] 2.500000e+02 3.220000e+02 4.070000e+02 6.120000e+02
## x_rep[7] 2.350000e+02 3.040000e+02 3.840000e+02 5.770000e+02
## x_rep[8] 2.250000e+02 2.900000e+02 3.660000e+02 5.510000e+02
## x_rep[9] 2.050000e+02 2.640000e+02 3.340000e+02 5.020000e+02
## x_rep[10] 1.850000e+02 2.390000e+02 3.020000e+02 4.540000e+02
## x_rep[11] 1.680000e+02 2.170000e+02 2.740000e+02 4.120000e+02
## x_rep[12] 1.600000e+02 2.070000e+02 2.610000e+02 3.940000e+02
## x_rep[13] 1.560000e+02 2.020000e+02 2.560000e+02 3.860000e+02
## x_rep[14] 1.560000e+02 2.020000e+02 2.550000e+02 3.840000e+02
## x_rep[15] 1.530000e+02 1.980000e+02 2.500000e+02 3.770000e+02
## x_rep[16] 1.530000e+02 1.980000e+02 2.500000e+02 3.770000e+02
## x_rep[17] 1.350000e+02 1.750000e+02 2.210000e+02 3.330000e+02
## x_rep[18] 1.290000e+02 1.670000e+02 2.110000e+02 3.190000e+02
## x_rep[19] 1.230000e+02 1.590000e+02 2.010000e+02 3.030000e+02
## x_rep[20] 1.100000e+02 1.430000e+02 1.810000e+02 2.740000e+02
## x_rep[21] 1.080000e+02 1.400000e+02 1.770000e+02 2.670000e+02
## x_rep[22] 1.030000e+02 1.340000e+02 1.700000e+02 2.570000e+02
## x_rep[23] 1.030000e+02 1.330000e+02 1.690000e+02 2.550000e+02
## x_rep[24] 9.800000e+01 1.280000e+02 1.620000e+02 2.440000e+02
## x_rep[25] 9.400000e+01 1.220000e+02 1.550000e+02 2.340000e+02
## x_rep[26] 7.900000e+01 1.030000e+02 1.310000e+02 1.980000e+02
## x_rep[27] 7.100000e+01 9.200000e+01 1.170000e+02 1.770000e+02
## x_rep[28] 3.500000e+01 4.600000e+01 5.900000e+01 9.000000e+01
## x_rep[29] 3.100000e+01 4.100000e+01 5.200000e+01 8.000000e+01
## x_rep[30] 3.100000e+01 4.100000e+01 5.200000e+01 8.000000e+01
## x_rep[31] 3.000000e+01 3.900000e+01 5.100000e+01 7.800000e+01
## x_rep[32] 2.100000e+01 2.800000e+01 3.700000e+01 5.700000e+01
## x_rep[33] 2.000000e+01 2.700000e+01 3.500000e+01 5.500000e+01
## x_rep[34] 1.800000e+01 2.400000e+01 3.100000e+01 4.800000e+01
## x_rep[35] 1.300000e+01 1.800000e+01 2.400000e+01 3.800000e+01
## x_rep[36] 1.200000e+01 1.700000e+01 2.200000e+01 3.500000e+01
## x_rep[37] 1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00
## lp__ -4.608346e+04 -4.608033e+04 -4.607747e+04 -4.607284e+04
## n_eff Rhat
## p[1] 251241.60 0.9999655
## p[2] 243459.13 1.0000150
## p[3] 263989.24 0.9999814
## p[4] 256446.81 0.9999685
## p[5] 258195.69 0.9999924
## p[6] 261859.80 0.9999669
## p[7] 261094.07 0.9999922
## p[8] 260180.84 0.9999726
## p[9] 264775.62 0.9999802
## p[10] 249604.24 0.9999645
## p[11] 257051.69 0.9999737
## p[12] 253736.18 0.9999704
## p[13] 266589.43 0.9999737
## p[14] 257040.97 0.9999665
## p[15] 259572.13 0.9999652
## p[16] 254715.59 0.9999632
## p[17] 251930.05 0.9999864
## p[18] 258953.12 0.9999597
## p[19] 256887.41 0.9999676
## p[20] 256434.43 0.9999876
## p[21] 252297.42 0.9999780
## p[22] 258416.77 0.9999771
## p[23] 242746.07 0.9999763
## p[24] 261879.31 0.9999661
## p[25] 259031.39 0.9999635
## p[26] 263775.84 0.9999669
## p[27] 252739.52 0.9999578
## p[28] 247500.75 0.9999896
## p[29] 257782.44 0.9999831
## p[30] 256642.67 0.9999830
## p[31] 258633.54 0.9999880
## p[32] 249934.73 0.9999619
## p[33] 256786.11 0.9999816
## p[34] 235431.19 0.9999968
## p[35] 235602.00 0.9999569
## p[36] 265217.57 0.9999708
## p[37] 241934.01 0.9999696
## phi 190092.57 1.0000202
## lambda 176396.68 1.0000104
## alpha 181940.70 0.9999950
## beta 176367.93 1.0000105
## x_rep[1] 145779.00 1.0000087
## x_rep[2] 145292.53 1.0000017
## x_rep[3] 145796.83 1.0000006
## x_rep[4] 145536.14 1.0000049
## x_rep[5] 145094.88 0.9999993
## x_rep[6] 145330.26 1.0000055
## x_rep[7] 144995.26 1.0000036
## x_rep[8] 145779.08 0.9999951
## x_rep[9] 145649.51 1.0000049
## x_rep[10] 145167.74 0.9999991
## x_rep[11] 145743.30 1.0000130
## x_rep[12] 145493.79 0.9999985
## x_rep[13] 145399.06 1.0000001
## x_rep[14] 145247.28 1.0000052
## x_rep[15] 145173.49 1.0000004
## x_rep[16] 145553.44 0.9999990
## x_rep[17] 145724.51 0.9999920
## x_rep[18] 145149.17 1.0000026
## x_rep[19] 145541.87 0.9999952
## x_rep[20] 145677.11 1.0000078
## x_rep[21] 144457.20 1.0000042
## x_rep[22] 144320.18 1.0000049
## x_rep[23] 145595.50 1.0000087
## x_rep[24] 145015.43 1.0000072
## x_rep[25] 145352.82 0.9999986
## x_rep[26] 145091.20 1.0000084
## x_rep[27] 145487.96 1.0000026
## x_rep[28] 144233.97 1.0000101
## x_rep[29] 144660.86 1.0000187
## x_rep[30] 145197.38 0.9999971
## x_rep[31] 145724.48 1.0000000
## x_rep[32] 144688.49 1.0000149
## x_rep[33] 143996.78 0.9999959
## x_rep[34] 143839.73 0.9999910
## x_rep[35] 145348.89 1.0000056
## x_rep[36] 145465.02 1.0000094
## x_rep[37] 142802.32 1.0000092
## lp__ 55502.94 1.0000325parameter estimates of beta distribution (for p)
phi1 = ss["phi","mean"] # posterior mean of phi
lambda1 = ss["lambda","mean"] # posterior mean of lambda
alpha1 = lambda1*phi1
beta1 = lambda1*(1-phi1)
beta_mean <- function(alpha,beta) alpha/(alpha+beta)
beta_var <- function(alpha,beta) alpha*beta/((alpha+beta)^2*(alpha+beta+1))phi1## [1] 0.00763126lambda1## [1] 1123.273alpha1## [1] 8.57199beta1## [1] 1114.701beta_mean(alpha1,beta1) # mean of p## [1] 0.00763126sqrt(beta_var(alpha1,beta1)) # sd of p## [1] 0.002595367pp = seq(0,1,length=1000)
plot(pp,dbeta(pp,alpha1,beta1),type="l") # distribution of p
extract posterior distributions (used later)
ex1 <- rstan::extract(fit1)
names(ex1)## [1] "p" "phi" "lambda" "alpha" "beta" "x_rep" "lp__"results
data points and posterior predictive distribution
p0=sum(data$x)/sum(data$n)
plot(data$n, data$x, type="n", xlab="投票総数", ylab="得票数", main="予測分布の95%区間", sub=paste("平均得票率=",signif(p0,4)),family=family)
#abline(a=0,b=phi1, col="red") # parameter estimates (=mean)
abline(a=0,b=p0, col="red") # = x/n
x_conf <- ss[paste("x_rep[",1:data$N,"]",sep=""), c("2.5%","97.5%","50%","mean")]
segments(data$n,x_conf[,1],data$n,x_conf[,2],col="green")
points(data$n, x_conf[,4],col="green" )
text(data$n, data$x, cex=0.5)
tail probability of posterior predictive distribution
ptail1 <- c()
for(i in 1:data$N) ptail1[i] = mean(ex1$x_rep[,i] <= data$x[i]) # lower tail
ptail1 <- 2*pmin(ptail1,1-ptail1) # posterior tail probability (two-tailed)
hist(ptail1,nclass=10,xlab="posterior probability")
plot of posterior tail probability
plot(data$n, ptail1, log="y", type="n", xlab="投票総数", ylab="事後確率", family=family)
text(data$n, ptail1);
abline(h=0.05,col="blue") # 5%
abline(h=0.01,col="orange") # 1%
abline(h=0.001,col="red") # 0.1%
VGAM includes beta-binomial distribution
library(VGAM)## Loading required package: stats4## Loading required package: splinestail probabilities of beta-binomial using the estimated parameter value
ptail2 = pbetabinom.ab(data$x,data$n,alpha1,beta1)
ptail2 <- 2*pmin(ptail2,1-ptail2) # p-value (two-tailed)
hist(ptail2,nclass=10,xlab="p-value") # check if uniform on (0,1)
plot(ptail1,ptail2); abline(0,1)
plot of frequntist p-value (based on the estimated parameters by stan)
plot(data$n, ptail2, log="y", type="n", xlab="投票総数", ylab="p値", family=family)
text(data$n, ptail2);
abline(h=0.05,col="blue") # 5%
abline(h=0.01,col="orange") # 1%
abline(h=0.001,col="red") # 0.1%
summary table of p-values
phi1## [1] 0.00763126lambda1## [1] 1123.273alpha1 ## [1] 8.57199beta1## [1] 1114.701cbind(as.character(data$city), ptail1, ptail2)## ptail1 ptail2
## [1,] "長岡市" "0.541486111111111" "0.533034497858158"
## [2,] "上越市" "0.000166666666666667" "1.35212852818374e-05"
## [3,] "中央区" "0.225236111111111" "0.211286126979409"
## [4,] "西区" "0.229458333333333" "0.216392919242591"
## [5,] "東区" "0.428916666666667" "0.419975328419225"
## [6,] "三条市" "0.413347222222222" "0.403675340388824"
## [7,] "新発田市" "0.667930555555555" "0.663417976154221"
## [8,] "柏崎市" "0.987916666666667" "0.989736273309983"
## [9,] "燕市" "0.524388888888889" "0.516590410785719"
## [10,] "秋葉区" "0.446888888888889" "0.438989658443659"
## [11,] "北区" "0.363944444444444" "0.355041559630298"
## [12,] "村上市" "0.811736111111111" "0.804060613268268"
## [13,] "南魚沼市" "0.829916666666667" "0.828793345096877"
## [14,] "十日町市" "0.686069444444444" "0.672559017970147"
## [15,] "佐渡市" "0.280694444444444" "0.262087718436888"
## [16,] "江南区" "0.644652777777778" "0.639326783608089"
## [17,] "西蒲区" "0.505597222222222" "0.497472845277565"
## [18,] "糸魚川市" "0.314013888888889" "0.295324349551976"
## [19,] "五泉市" "0.817555555555556" "0.816358344969879"
## [20,] "阿賀野市" "0.861319444444445" "0.861162249253763"
## [21,] "魚沼市" "0.784805555555556" "0.775499452507736"
## [22,] "見附市" "0.412902777777778" "0.402690720144685"
## [23,] "小千谷市" "0.950125" "0.952145741735418"
## [24,] "南区" "0.338263888888889" "0.328598258984664"
## [25,] "妙高市" "0.931875" "0.925371162791766"
## [26,] "胎内市" "0.946361111111111" "0.949823252158155"
## [27,] "加茂市" "0.767666666666667" "0.766145139672933"
## [28,] "阿賀町" "0.331847222222222" "0.314430311109931"
## [29,] "聖籠町" "0.689097222222222" "0.67871724054185"
## [30,] "田上町" "0.543152777777778" "0.527886110419607"
## [31,] "津南町" "0.898708333333333" "0.900298486012659"
## [32,] "弥彦村" "0.721708333333333" "0.721820651192607"
## [33,] "湯沢町" "0.598291666666667" "0.592661502755018"
## [34,] "関川村" "0.175638888888889" "0.161352769588204"
## [35,] "刈羽村" "0.336486111111111" "0.321867690828018"
## [36,] "出雲崎町" "0.8475" "0.839244845842987"
## [37,] "粟島浦村" "0.987805555555556" "0.989887001551361"fitting the ordinary binomial distribution (THIS IS NOT A VALID ANALYSIS)
p
(x = sum(data$x))## [1] 8009(n = sum(data$n))## [1] 1004439(p= x/n)## [1] 0.007973605tail probabilities of observed values (p-values) : DONT USE THIS!!!
ptail3 = pbinom(data$x,data$n,p)
ptail3 <- 2*pmin(ptail3,1-ptail3) # p-value (two-tailed) with Binomial
hist(ptail3,nclass=10,xlab="p-value (Binomial)") # check if uniform on (0,1)
plot(ptail2,ptail3); abline(0,1)
ptail3## [1] 4.928394e-05 1.205873e-190 0.000000e+00 2.220446e-16 3.494165e-05
## [6] 8.329335e-05 2.177306e-01 1.745646e-01 1.492489e-02 2.293801e-03
## [11] 1.707704e-04 1.464160e-02 9.129030e-01 1.167704e-03 2.254622e-10
## [16] 2.373768e-01 3.139326e-02 2.846726e-08 9.846748e-01 7.925077e-01
## [21] 3.008091e-02 9.151586e-03 4.668626e-01 1.572925e-03 1.944645e-01
## [26] 5.306942e-01 8.182928e-01 2.954428e-03 1.313150e-01 4.561573e-02
## [31] 8.169958e-01 7.616075e-01 4.940411e-01 3.254725e-03 4.256557e-02
## [36] 4.870894e-01 9.225547e-01posterior predictive distributions of all x’s
for(i in 1:data$N) {
hist(ex1$x_rep[,i], nclass=100, xlab="得票数", main=paste("事後予測分布 (",data$city[i],")"), sub=paste(i,":", data$city[i],", 事後確率=", signif(ptail1[i],4), ", p値=", signif(ptail2[i],4)), family= family) # histogram of posterior distribution
abline(v = data$x[i], col="red" ) # observed value
}
