# EXAMPLE OF BAYESIAN INFERENCE USING TRAPEZOIDAL 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, nsum, 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. The likelihood is # indeterminate for the model with zeros when lambda is zero; the value # returned is arbitrarily set to zero. loglike.with.zeros <- function(lambda,nsum,k) { nsum * log(lambda) - k * lambda } loglike.no.zeros <- function(lambda,nsum,k) { if (lambda==0) { return (0) } nsum * 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,nsum,k) { exp (loglike.with.zeros(lambda,nsum,k) - loglike.with.zeros(nsum/k,nsum,k)) } upost.no.zeros <- function (lambda,nsum,k) { exp (loglike.no.zeros(lambda,nsum,k) - loglike.with.zeros(nsum/k,nsum,k)) } # COMPUTE VARIOUS THINGS FOR THE TWO MODELS FOR A GIVEN DATA SET. do.it <- function (n) { cat("The data points:",n,"\n\n") k <- length(n) nsum <- sum(n) cat("Number of points:",k,"\n") cat("Sum of points:",nsum,"\n") cat("Mean of points:",nsum/k,"\n\n") norm.with.zeros <- trapezoidal.rule (upost.with.zeros, 0, 3, 10, nsum, k) norm.no.zeros <- trapezoidal.rule (upost.no.zeros, 0, 3, 10, nsum, 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 <- trapezoidal.rule ( function (lambda,nsum,k) lambda * upost.with.zeros(lambda,nsum,k), 0, 3, 10, nsum, k) / norm.with.zeros E.lambda.no.zeros <- trapezoidal.rule ( function (lambda,nsum,k) lambda * upost.no.zeros(lambda,nsum,k), 0, 3, 10, nsum, 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 TRAPEZOID RULE. The arguments are the # function to integrate, the low and high bounds, and the number of times # to double the number of points (so the final number of points used is # 2^K+1). Additional arguments to pass to the function being integrated # can follow. Results are printed for each approximations along the way, # so you can see how fast the integral converges. trapezoidal.rule <- function(f,a,b,K,...) { # Initial estimate based on just the two end points. h <- b-a n <- 1 I <- (f(a,...)+f(b,...))*(h/2) cat("Trapezoidal integration:\n") cat(0,I,"\n") for (j in 1:K) { t <- a+(h/2) p <- 0 for (i in 0:(n-1)) { p <- p + f(t+i*h,...) } I <- (I+h*p)/2 cat(j,I,"\n") h <- h/2 n <- 2*n } cat("\n") I }