# DEMO OF SIMPLE MONTE CARLO INFERENCE FOR A BAYESIAN MODEL.
#
# Radford M. Neal, January 2013.
#
# We observe a pair of non-negative integer data points, y1 and y2,
# which we model as being drawn independently from binomial
# distributions with the same n (a known constant) and probability
# parameter theta1 for y1 and theta2 for y2.
#
# Two prior distributions for (theta1,theta2) are considered. In the
# "unordered" prior, theta1 and theta2 are independent, both with the
# the uniform(0,1) distribution - ie, (theta1,theta2) is uniform over
# the square where theta1 and theta2 range from 0 to 1. In the
# "ordered" prior, (theta1,theta2) is uniform over the triangle in
# which theta1 and theta2 range from 0 to 1 with theta2 >= theta1.
# SAMPLE FROM THE PRIOR DISTRIBUTION. The argument 'N' is the number
# of (theta1,theta2) values to sample from the prior. The argument
# 'ordered' is TRUE for the ordered prior and FALSE for the unordered
# prior (default is FALSE). The value is a matrix with two columns
# and 'N' rows, with each row being an independent sample from the prior.
sample_prior <- function (N, ordered=FALSE)
{
theta <- matrix (NA, N, 2)
for (i in 1:N)
{ th <- runif(2)
if (ordered) theta[i,] <- sort(th)
else theta[i,] <- th
}
theta
}
# SAMPLE FROM THE POSTERIOR DISTRIBUTION. This is done by sampling
# from the prior for (theta1,theta2), then sampling (y1,y2) using
# these values for theta1 and theta2, and finally discarding those
# samples that don't match the observed y1 and y2.
#
# The argument 'N' is the number of values to sample from the prior.
# The argument 'n' is the known n parameter of the binomial
# distributions. Arguments 'y1' and 'y2' are the observed data. Any
# further arguments are passed on to 'sample_prior'.
#
# The value is a matrix with two columns and some number of rows (from
# 0 to 'N'), containing (theta1,theta2) values drawn independently
# from the posterior distribution. An estimate of the marginal
# likelihood for the model (which is just the fraction of prior values
# that match the observed data) is printed.
sample_posterior <- function (N, n, y1, y2, ...)
{
theta <- sample_prior(N,...)
prior_y1 <- rbinom(N,n,theta[,1])
prior_y2 <- rbinom(N,n,theta[,2])
matches <- prior_y1==y1 & prior_y2==y2
cat("Marginal likelihood:",mean(matches),"\n")
theta[matches,,drop=FALSE]
}
# PLOT PAIRS OF THETA VALUES. Produces a scatterplot of (theta1,theta2)
# values, with a diagonal line drawn where theta1=theta2. The argument
# is a matrix with two columns containing values for theta1 and theta2.
plot_thetas <- function (theta)
{
plot (theta[,1], theta[,2], xlim=c(0,1), ylim=c(0,1), pch=20)
abline(0,1)
}