# Solution to STA 410/2102 Assignment #4, Spring 2003 - Program code.


# LOG LIKELIHOOD FUNCTION.  The data is given in terms of its sufficiet 
# statistics, the sum of data points, dsum, and the number of points, n.
# The factorial factors in the likelihood are omitted.

log.likelihood <- function(theta,dsum,n) 
{ dsum * theta - n * exp(theta)
}


# COMPUTE THE UNNORMALIZED POSTERIOR DENSITY.  The log likelihood is offset by 
# the log likelihood for the model with the sample mean as the Poisson mean, 
# to prevent underflow problems (except that this is not done if the sample
# mean is zero).  The parameters are the value for the parameter theta, the
# sum of the data points, and the number of data points.

uposterior <- function (theta,dsum,n) 
{ 
  if (dsum>0)
  { offset <- log.likelihood(log(dsum/n),dsum,n)
  }
  else
  { offset <- 0
  }

  return (dnorm(theta,0.7,1.5) * exp (log.likelihood(theta,dsum,n) - offset))
}


# COMPUTE VARIOUS THINGS FOR A GIVEN DATA SET.

do.it <- function (d, N=100)
{ 

  cat("\nNumber of function evaluations to use for integration:",N,"\n\n")

  cat("The data points:",d,"\n\n")

  n <- length(d)
  dsum <- sum(d)

  cat("Number of data points:",n,"\n")
  cat("Sum of data points:",dsum,"\n")
  cat("Mean of data points:",dsum/n,"\n")
  cat("Probability of a zero using mean for lambda:",exp(-dsum/n),"\n\n");

  norm <- inf.midpoint (uposterior, -Inf, +Inf, N, dsum, n)

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

  E.theta <- (1/norm) * inf.midpoint (
     function (theta,dsum,k) theta * uposterior(theta,dsum,n),
     -Inf, +Inf, N, dsum, n)

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

  P.pos <- (1/norm) * inf.midpoint (uposterior, 0, +Inf, N, dsum, n)

  cat("Probability that theta is positive:", P.pos, "\n");

  pred0 <- (1/norm) * inf.midpoint (
     function (theta,dsum,k) exp(-exp(theta)) * uposterior(theta,dsum,n),
     -Inf, +Inf, N, dsum, n)

  cat("Predictive probability of a zero:", pred0, "\n");
}


# NUMERICAL INTEGRATION FROM (POSSIBLY) -INFINITY TO +INFINITY.  The arguments
# are the function to integrate, the bounds of integration, which may be 
# -Inf or +Inf, the number of points to evaluate the function at, and any
# additional arguments to the function.  The integration is done by transforming
# from (-Inf,+Inf) to (0,1) by means of 1/(1+exp(-x)).

inf.midpoint <- function (f,a,b,N,...)
{ ta <- 1/(1+exp(-a))
  tb <- 1/(1+exp(-b))
  midpoint.rule (function (y,...) f(log(y/(1-y)),...) / (y*(1-y)), 
                 ta, tb, 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
}