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)