Revised binomial examples

johnmyleswhite · Mar 22, 2012 · 23984d3 · 23984d3
1 parent 2caedf5
commit 23984d3
Show file tree

Hide file tree

Showing 5 changed files with 43 additions and 2 deletions.
diff --git a/graphs/binomial/broken_plot.png b/graphs/binomial/broken_plot.png
diff --git a/graphs/binomial/density.png b/graphs/binomial/density.png
diff --git a/graphs/binomial/plot.png b/graphs/binomial/plot.png
diff --git a/scripts/binomial/binomial.R b/scripts/binomial/binomial.R
@@ -1,8 +1,13 @@
+# Basic inference.
 library('rjags')
 
-df <- read.csv(file.path('data', 'binomial', 'binomial.csv'))
+df <- read.csv(file.path('data',
+                         'binomial',
+                         'binomial.csv'))
 
-jags <- jags.model(file.path('bugs', 'binomial', 'binomial.bugs'),
+jags <- jags.model(file.path('bugs',
+                             'binomial',
+                             'binomial.bugs'),
                    data = list('x' = with(df, X),
                                'N' = nrow(df)),
                    n.chains = 4,
@@ -12,10 +17,13 @@ mcmc.samples <- coda.samples(jags,
                              c('p'),
                              2500)
 
+png(file.path('graphs', 'binomial', 'plot.png'))
 plot(mcmc.samples)
+dev.off()
 
 summary(mcmc.samples)
 
+# Compare approximate posterior mean with analytic posterior mean.
 alpha <- with(df, sum(X) + 1)
 beta <- nrow(df) + 2 - alpha
 
@@ -30,3 +38,14 @@ 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'))
+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
+credible.interval <- c(qbeta(0.025, alpha, beta),
+                       qbeta(0.975, alpha, beta))
+credible.interval
+binom.test(with(df, sum(X)) + 1, nrow(df) + 2)$conf.int
diff --git a/scripts/binomial/broken_binomial.R b/scripts/binomial/broken_binomial.R
@@ -0,0 +1,22 @@
+# Basic inference with no adaptive phase and few posterior samples.
+library('rjags')
+
+df <- read.csv(file.path('data',
+                         'binomial',
+                         'binomial.csv'))
+
+jags <- jags.model(file.path('bugs',
+                             'binomial',
+                             'binomial.bugs'),
+                   data = list('x' = with(df, X),
+                               'N' = nrow(df)),
+                   n.chains = 4,
+                   n.adapt = 0)
+
+mcmc.samples <- coda.samples(jags,
+                             c('p'),
+                             100)
+
+png(file.path('graphs', 'binomial', 'broken_plot.png'))
+plot(mcmc.samples)
+dev.off()