# STA 410, ASSIGNMENT 4, SPRING 2004.


# INTEGRATE A ONE-DIMENSIONAL FUNCTION USING SIMPSON'S RULE.  Arguments
# are the function (which may take extra arguments beyond the one being
# integrated over), the lower limit of integration, the upper limit of
# integration, the number of intervals for the approximation (the number
# of points used is twice this plus one), and any extra arguments to
# pass on to the function.  The value is the approximation to the integral.

simpson <- function (f, a, b, n, ...)
{
  if (n<1 || a>b) 
  { stop("Invalid argument")
  }

  h = (b-a)/(2*n)

  s <- f(a,...)
  for (i in 1:n)
  { s <- s + 4*f(a+(2*i-1)*h)
  }
  if (n>1)
  { for (i in 1:(n-1))
    { s <- s + 2*f(a+2*i*h)
    }
  }
  s <- s + f(b,...)

  s * h/3
}


# LIKELIHOOD FUNCTION.  Takes the x and y data vectors as arguments, along
# with values for alpha and beta, and returns the probability of y given
# alpha, beta, and x.

likelihood <- function (x, y, alpha, beta)
{
  prod (dpois(y,alpha+beta*x))
}


# COMPUTE MARGINAL LIKELIHOOD AND E(ALPHA) FOR H1.  Takes the x and y
# data vectors as arguments, along with the number of number of intervals
# to use for Simpson's Rule.  Returns a list whose elements are the marginal
# likelihood (marg.lik) and the expectation (E1).

H1 <- function (x, y, n)
{
  marg.lik <- 
    simpson (function (alpha) likelihood(x,y,alpha,0) / 10, 0, 10, n)
  E1 <- (1/marg.lik) *
    simpson (function (alpha) alpha * likelihood(x,y,alpha,0) / 10, 0, 10, n)

  list (marg.lik=marg.lik, E1=E1)
}


# COMPUTE MARGINAL LIKELIHOOD AND E(BETA) FOR H2.  Takes the x and y
# data vectors as arguments, along with the number of number of intervals
# to use for Simpson's Rule.  Returns a list whose elements are the marginal
# likelihood (marg.lik) and the expectation (E1).

H2 <- function (x, y, n)
{
  marg.lik <- 
    simpson (function (beta) likelihood(x,y,0,beta) / 5, 0, 5, n)
  E1 <- (1/marg.lik) *
    simpson (function (beta) beta * likelihood(x,y,0,beta) / 5, 0, 5, n)

  list (marg.lik=marg.lik, E1=E1)
}


# COMPUTE MARGINAL LIKELIHOOD AND E(ALPHA+BETA) FOR H3.  Takes the x and y
# data vectors as arguments, along with the number of number of intervals
# to use for Simpson's Rule.  Returns a list whose elements are the marginal
# likelihood (marg.lik) and the expectation (E1).

H3 <- function (x, y, n)
{
  marg.lik <- 
   simpson (function (alpha) simpson 
              (function (beta) likelihood(x,y,alpha,beta) / (5*4), 0, 5, n),
            0, 4, n)
             
  E1 <- (1/marg.lik) *
   simpson (function (alpha) simpson 
              (function (beta) (alpha+beta) * likelihood(x,y,alpha,beta) / (5*4), 
               0, 5, n),
            0, 4, n)

  list (marg.lik=marg.lik, E1=E1)
}