binomial

Basic inference.

library('rjags')

## Loading required package: coda
## Linked to JAGS 3.4.0
## Loaded modules: basemod,bugs

library(ggplot2)

df <- data.frame(X=c(0,0,0,1,0))

binomial_code <- '
  model
    {
    for (i in 1:N)
    {
    x[i] ~ dbern(p)
    }
    
    p ~ dbeta(alpha, beta)

    alpha <- 1
    beta <- 1
    }
'

writeLines(binomial_code,con="binomial.txt")

jags <- jags.model("binomial.txt",
                   data = list('x' = with(df, X),
                               'N' = nrow(df)),
                   n.chains = 4,
                   n.adapt = 1000)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 9
## 
## Initializing model

mcmc.samples <- coda.samples(jags,
                             c('p'),
                             5000)

#png(file.path('graphs', 'binomial', 'plot.png'))
plot(mcmc.samples)

dev.off()

## null device 
##           1

summary(mcmc.samples)

## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 4 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##       0.284957       0.160301       0.001133       0.001142 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.0439 0.1587 0.2638 0.3898 0.6393

# Compare approximate posterior mean with analytic posterior mean.
alpha <- with(df, sum(X) + 1)
beta <- nrow(df) + 2 - alpha

alpha / (alpha + beta)

## [1] 0.2857143

# Compare distribution of samples with theoretical density.
samples <- data.frame(MCMC = as.matrix(mcmc.samples)[, 1])
density.plot <- ggplot(samples, aes(x = MCMC)) +
  geom_density(aes(color = 'Sampling'))
density.plot 

x <- seq(0, 1, by = 0.01)
analytic <- data.frame(x = x, y = dbeta(x, alpha, beta))
density.plot <- density.plot +
  geom_line(data = analytic, aes(x = x, y = y, color = 'Analytic'))
density.plot

#ggsave(file.path('graphs', 'binomial', 'density.png'))

# Compare confidence intervals with credible intervals.
# binom.test(1, 5)
# binom.test(2, 7)

binom.test(with(df, sum(X)), nrow(df))$conf.int

## [1] 0.005050763 0.716417936
## attr(,"conf.level")
## [1] 0.95

credible.interval <- c(qbeta(0.025, alpha, beta),
                       qbeta(0.975, alpha, beta))
credible.interval

## [1] 0.04327187 0.64123458

binom.test(with(df, sum(X)) + 1, nrow(df) + 2)$conf.int

## [1] 0.03669257 0.70957914
## attr(,"conf.level")
## [1] 0.95