# EXAMPLE OF BAYESIAN INFERENCE USING MIDPOINT RULE INTEGRATION. # # Bayesian inference for a Poisson(lambda) model with prior for lambda uniform # over [0,3]. An alternate model says that zeros eliminated before the data # was seen. # LOG LIKELIHOOD FUNCTIONS FOR THE TWO MODELS. The data is given in terms # of its sufficiet statistics, the sum of data points, dsum, and the number # of points, k. The factorial factors in the likelihood are omitted, but this # is OK since they are the same for the two models. Note that the likelihood # is undefined for the model with no zeros when lambda is zero. loglike.with.zeros <- function(lambda,dsum,k) { dsum * log(lambda) - k * lambda } loglike.no.zeros <- function(lambda,dsum,k) { dsum * log(lambda) - k * (lambda + log(1-exp(-lambda))) } # COMPUTE THE UNNORMALIZED POSTERIOR DENSITIES. The likelihood is offset by # the likelihood for the model with zeros for with the sample mean as the # parameter, to prevent underflow problems. Since the prior is uniform, it # is omitted here. upost.with.zeros <- function (lambda,dsum,k) { exp (loglike.with.zeros(lambda,dsum,k) - loglike.with.zeros(dsum/k,dsum,k)) } upost.no.zeros <- function (lambda,dsum,k) { exp (loglike.no.zeros(lambda,dsum,k) - loglike.with.zeros(dsum/k,dsum,k)) } # COMPUTE VARIOUS THINGS FOR THE TWO MODELS FOR A GIVEN DATA SET. do.it <- function (d, n=1000) { cat("The data points:",d,"\n\n") k <- length(d) dsum <- sum(d) cat("Number of points:",k,"\n") cat("Sum of points:",dsum,"\n") cat("Mean of points:",dsum/k,"\n\n") norm.with.zeros <- midpoint.rule (upost.with.zeros, 0, 3, n, dsum, k) norm.no.zeros <- midpoint.rule (upost.no.zeros, 0, 3, n, dsum, k) cat("Normalizing constants for posterior:", norm.with.zeros,"with zeros,", norm.no.zeros,"no zeros\n") cat("Bayes factor (with zeros / no zeros):", norm.with.zeros/norm.no.zeros,"\n\n") E.lambda.with.zeros <- midpoint.rule ( function (lambda,dsum,k) lambda * upost.with.zeros(lambda,dsum,k), 0, 3, n, dsum, k) / norm.with.zeros E.lambda.no.zeros <- midpoint.rule ( function (lambda,dsum,k) lambda * upost.no.zeros(lambda,dsum,k), 0, 3, n, dsum, k) / norm.no.zeros cat("Posterior means for lambda:", E.lambda.with.zeros,"with zeros,", E.lambda.no.zeros,"no zeros\n") } # NUMERICAL INTEGRATION USING THE MIDPOINT RULE. The arguments are the # function to integrate, the low and high bounds, and the number of points # to use. Additional arguments to pass to the function being integrated # can follow. midpoint.rule <- function(f,a,b,n,...) { h <- (b-a)/n s <- 0 for (i in 1:n) { s <- s+f(a+h*(i-0.5),...) } h*s }