20190726nigata-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))remove jouetsu from data : THIS IS DIFFERENT FROM THE PREVIOUS ANALYSIS
In data2 (N, n, x), joetsu is removed from data. The original data is now renamed as N_new, n_new, x_new. N, n, x are used for computing posterior distribution. Then predictive distribution is computed for n_new, and posterior probability is computed at x_new.
id2 = (1:data$N)[-2] # 2=joetsu
data2 = list(N=data$N-1, x=data$x[id2], n=data$n[id2], id=id2, city=data$city[id2], N_new=data$N, x_new=data$x, n_new=data$n, city_new=data$city)plot
plot(data2$n_new, data2$x_new, type="n", xlab="投票総数", ylab="得票数", family=family); text(data2$n_new, data2$x_new, cex=0.5,col="red") ; text(data2$n, data2$x, data2$id, 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
betabinom3 = stan_model("betabinom3.stan")
fit1 <- sampling(betabinom3, data=data2, chains=8, warmup=2000, iter=20000)stan diagnostics
fit1## Inference for Stan model: betabinom3.
## 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.01 0.00 0.00 0.01 0.01 0.01 0.01
## 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.00 0.00 0.01 0.01
## p[15] 0.01 0.00 0.00 0.01 0.01 0.01 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.00 0.00 0.01 0.01
## p[18] 0.01 0.00 0.00 0.01 0.01 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.00 0.01 0.01 0.01
## p[28] 0.01 0.00 0.00 0.01 0.01 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.01 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.00 0.00 0.01 0.01
## p[34] 0.01 0.00 0.00 0.00 0.01 0.01 0.01
## p[35] 0.01 0.00 0.00 0.01 0.01 0.01 0.01
## p[36] 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 2975.21 2.77 879.03 1573.74 2347.95 2868.29 3482.28
## alpha 23.58 0.02 7.03 12.42 18.56 22.70 27.61
## beta 2951.63 2.75 872.07 1561.07 2329.48 2845.77 3454.83
## x_rep[1] 934.65 0.54 205.69 572.00 793.00 921.00 1060.00
## x_rep[2] 684.96 0.40 151.32 418.00 581.00 674.00 777.00
## x_rep[3] 597.47 0.35 132.34 364.00 506.00 588.00 678.00
## x_rep[4] 525.36 0.31 116.62 319.00 445.00 517.00 597.00
## x_rep[5] 429.31 0.25 95.68 260.00 363.00 423.00 488.00
## x_rep[6] 348.40 0.21 78.09 211.00 294.00 343.00 396.00
## x_rep[7] 328.85 0.19 73.71 199.00 278.00 324.00 374.00
## x_rep[8] 313.24 0.18 70.44 189.00 265.00 308.00 356.00
## x_rep[9] 285.63 0.17 64.39 172.00 241.00 281.00 325.00
## x_rep[10] 258.42 0.15 58.38 155.00 218.00 254.00 294.00
## x_rep[11] 234.28 0.14 53.22 140.00 197.00 231.00 267.00
## x_rep[12] 223.64 0.13 50.89 134.00 189.00 220.00 255.00
## x_rep[13] 218.83 0.13 49.78 131.00 184.00 215.00 249.00
## x_rep[14] 218.38 0.13 49.76 130.00 184.00 215.00 249.00
## x_rep[15] 214.09 0.13 48.83 128.00 180.00 211.00 244.00
## x_rep[16] 213.71 0.13 48.73 128.00 180.00 210.00 244.00
## x_rep[17] 189.07 0.11 43.34 113.00 159.00 186.00 216.00
## x_rep[18] 180.81 0.11 41.55 107.00 152.00 178.00 206.00
## x_rep[19] 171.81 0.10 39.59 102.00 144.00 169.00 196.00
## x_rep[20] 154.86 0.09 35.95 91.00 130.00 152.00 177.00
## x_rep[21] 151.56 0.09 35.21 90.00 127.00 149.00 173.00
## x_rep[22] 145.21 0.09 33.82 86.00 122.00 143.00 166.00
## x_rep[23] 144.32 0.09 33.61 85.00 121.00 142.00 165.00
## x_rep[24] 138.00 0.08 32.26 81.00 116.00 136.00 158.00
## x_rep[25] 132.42 0.08 31.00 78.00 111.00 130.00 151.00
## x_rep[26] 111.47 0.07 26.40 65.00 93.00 110.00 128.00
## x_rep[27] 99.69 0.06 23.86 58.00 83.00 98.00 114.00
## x_rep[28] 49.76 0.03 12.90 27.00 41.00 49.00 58.00
## x_rep[29] 43.99 0.03 11.63 24.00 36.00 43.00 51.00
## x_rep[30] 43.93 0.03 11.59 24.00 36.00 43.00 51.00
## x_rep[31] 42.88 0.03 11.39 23.00 35.00 42.00 50.00
## x_rep[32] 30.76 0.02 8.67 16.00 25.00 30.00 36.00
## x_rep[33] 29.68 0.02 8.46 15.00 24.00 29.00 35.00
## x_rep[34] 25.93 0.02 7.61 13.00 21.00 25.00 31.00
## x_rep[35] 19.83 0.02 6.19 9.00 15.00 19.00 24.00
## x_rep[36] 18.43 0.02 5.87 8.00 14.00 18.00 22.00
## x_rep[37] 1.79 0.00 1.39 0.00 1.00 2.00 3.00
## lp__ -45408.98 0.02 4.53 -45418.77 -45411.83 -45408.64 -45405.76
## 97.5% n_eff Rhat
## p[1] 0.01 201705 1
## p[2] 0.01 187530 1
## p[3] 0.01 197439 1
## p[4] 0.01 200015 1
## p[5] 0.01 205708 1
## p[6] 0.01 204975 1
## p[7] 0.01 198228 1
## p[8] 0.01 193472 1
## p[9] 0.01 197279 1
## p[10] 0.01 193216 1
## p[11] 0.01 198090 1
## p[12] 0.01 203409 1
## p[13] 0.01 187908 1
## p[14] 0.01 187676 1
## p[15] 0.01 199284 1
## p[16] 0.01 196178 1
## p[17] 0.01 174493 1
## p[18] 0.01 185510 1
## p[19] 0.01 197305 1
## p[20] 0.01 196491 1
## p[21] 0.01 201517 1
## p[22] 0.01 193835 1
## p[23] 0.01 197673 1
## p[24] 0.01 206973 1
## p[25] 0.01 204791 1
## p[26] 0.01 223610 1
## p[27] 0.01 173517 1
## p[28] 0.01 177553 1
## p[29] 0.01 185722 1
## p[30] 0.01 189775 1
## p[31] 0.01 205773 1
## p[32] 0.01 197978 1
## p[33] 0.01 162531 1
## p[34] 0.01 161178 1
## p[35] 0.01 192443 1
## p[36] 0.01 188522 1
## phi 0.01 140068 1
## lambda 5002.43 100365 1
## alpha 39.78 97174 1
## beta 4962.91 100397 1
## x_rep[1] 1379.00 145301 1
## x_rep[2] 1014.00 145101 1
## x_rep[3] 885.00 145363 1
## x_rep[4] 779.00 145845 1
## x_rep[5] 636.00 145183 1
## x_rep[6] 518.00 145001 1
## x_rep[7] 489.00 144927 1
## x_rep[8] 467.00 145737 1
## x_rep[9] 425.00 145420 1
## x_rep[10] 385.00 145620 1
## x_rep[11] 350.00 145306 1
## x_rep[12] 334.00 145574 1
## x_rep[13] 326.00 145668 1
## x_rep[14] 326.00 145406 1
## x_rep[15] 320.00 145210 1
## x_rep[16] 319.00 145012 1
## x_rep[17] 283.00 145327 1
## x_rep[18] 271.00 145111 1
## x_rep[19] 257.00 145561 1
## x_rep[20] 233.00 144994 1
## x_rep[21] 228.00 145832 1
## x_rep[22] 219.00 145311 1
## x_rep[23] 217.00 145422 1
## x_rep[24] 208.00 145234 1
## x_rep[25] 200.00 144683 1
## x_rep[26] 169.00 145394 1
## x_rep[27] 152.00 145942 1
## x_rep[28] 78.00 143969 1
## x_rep[29] 69.00 144857 1
## x_rep[30] 69.00 145239 1
## x_rep[31] 68.00 145570 1
## x_rep[32] 49.00 142616 1
## x_rep[33] 48.00 143532 1
## x_rep[34] 42.00 144524 1
## x_rep[35] 33.00 145362 1
## x_rep[36] 31.00 143997 1
## x_rep[37] 5.00 142930 1
## lp__ -45401.10 57086 1
##
## Samples were drawn using NUTS(diag_e) at Fri Jul 26 11:29:36 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)