# EXAMPLE OF BAYESIAN INFERENCE USING MIDPOINT RULE INTEGRATION.
#
# Bayesian inference for a Poisson(lambda) model with prior for lambda uniform 
# over [0,3]. 


# LOG LIKELIHOOD FUNCTION.  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.

log.likelihood <- function(lambda,dsum,k) 
{ dsum * log(lambda) - k * lambda 
}


# COMPUTE THE UNNORMALIZED POSTERIOR DENSITY.  The log likelihood is offset by 
# the log likelihood for the model with the sample mean as the parameter, to 
# prevent underflow problems.  Since the prior is uniform, it can be omitted,
# except that zero is returned outside the range [0,3].

uposterior <- function (lambda,dsum,k) 
{ if (lambda<0 || lambda>3)
  { return (0)
  }
  else
  { return (exp (log.likelihood(lambda,dsum,k) - log.likelihood(dsum/k,dsum,k)))
  }
}


# COMPUTE VARIOUS THINGS FOR FOR A GIVEN DATA SET.

do.it <- function (d, n=100)
{ 
  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")
  cat("Probability of a zero using mean for lambda:",exp(-dsum/k),"\n\n");

  norm <- midpoint.rule (uposterior, 0, 3, n, dsum, k)

  cat("Normalizing constant for posterior:",norm,"\n");

  E.lambda <- (1/norm) * midpoint.rule (
     function (lambda,dsum,k) lambda * uposterior(lambda,dsum,k),
     0, 3, n, dsum, k)

  cat("Posterior mean for lambda:", E.lambda, "\n");

  pred0 <- (1/norm) * midpoint.rule (
     function (lambda,dsum,k) exp(-lambda) * uposterior(lambda,dsum,k),
     0, 3, n, dsum, k)

  cat("Predictive probability of a zero:", pred0, "\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
}