Computational Actuarial Science - Ex15

library(R2WinBUGS)
library(lattice)

hospital <- data.frame(n=c(47,148,119,810,211,196,148,215,207,97,256),r=c(0,18,8,46,8,13,9,31,14,8,29))

# preparing everything

k <- nrow(hospital)
n <- hospital$n
r <- hospital$r
data1 <- list("k", "n", "r") #List of names

inits <- function()
{
 list(p = rep(0.1,11))
}

# running the MCMC model

sim<-bugs(data1,inits,model.file = "inst ranking 1.txt",parameters=c("p")
          ,n.chains = 1,n.iter= 1000, bugs.directory = "C:/Users/natan/Google Drive/Act/Computational Actuarial Science/ HW/15/WinBUGS14")



post<-as.mcmc.list(sim)
summary(post)

## 
## Iterations = 501:1000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##              Mean       SD  Naive SE Time-series SE
## deviance 56.70422 5.045022 0.2256202      0.2256202
## p[1]      0.01984 0.018859 0.0008434      0.0008434
## p[10]     0.09283 0.029864 0.0013356      0.0013356
## p[11]     0.11752 0.020219 0.0009042      0.0008549
## p[2]      0.12806 0.027596 0.0012341      0.0012341
## p[3]      0.07184 0.023707 0.0010602      0.0010602
## p[4]      0.05797 0.008458 0.0003782      0.0003782
## p[5]      0.04201 0.013110 0.0005863      0.0005863
## p[6]      0.07091 0.017889 0.0008000      0.0007334
## p[7]      0.06699 0.019649 0.0008787      0.0008787
## p[8]      0.14575 0.024640 0.0011019      0.0011019
## p[9]      0.07243 0.017616 0.0007878      0.0007878
## 
## 2. Quantiles for each variable:
## 
##               2.5%      25%      50%      75%    97.5%
## deviance 4.901e+01 52.88500 56.04500 59.70000 68.31975
## p[1]     7.821e-04  0.00620  0.01303  0.02839  0.06576
## p[10]    4.443e-02  0.07209  0.08901  0.11060  0.16072
## p[11]    8.304e-02  0.10410  0.11600  0.12843  0.16256
## p[2]     8.138e-02  0.10837  0.12535  0.14402  0.18632
## p[3]     3.456e-02  0.05529  0.06910  0.08516  0.12910
## p[4]     4.279e-02  0.05204  0.05697  0.06383  0.07466
## p[5]     1.995e-02  0.03238  0.04101  0.05020  0.06914
## p[6]     3.889e-02  0.05795  0.06989  0.08226  0.10906
## p[7]     3.328e-02  0.05361  0.06531  0.07922  0.11149
## p[8]     1.012e-01  0.12950  0.14380  0.16173  0.19943
## p[9]     4.347e-02  0.06007  0.07120  0.08197  0.11054

plot(post)

xyplot(post)

densityplot(post)

acfplot(post)

cumuplot(post)