# STA 410/2102, Assignment #3.  R. M. Neal.


# CONSTANTS PERTAINING TO THE MODEL.

low <- 0
high <- 10
sigma <- 1


# THE LIKELIHOOD FUNCTION.
#
# This is done in a simple fashion, but it would be more efficient to
# pre-compute the sample means and use them as sufficient statistics.

likelihood <- function (theta1,theta2,data1,data2)
{ 
  prod(exp(-(data1-theta1)^2/(2*sigma^2))) *
    prod(exp(-(data2-theta2)^2/(2*sigma^2)))
}


# FIND THE POSTERIOR MEANS OF THE PARAMETERS.  Passed the data vectors and
# the number of points to use in the trapezoidal rule.  Returns the vector
# of the two posterior means.

posterior.means <- function (data1,data2,n)
{
  norm <- triangular.int(likelihood,low,high,n,data1,data2)

  mean1 <- triangular.int(
    function(theta1,theta2,data1,data2) 
      theta1 * likelihood(theta1,theta2,data1,data2),
    low, high, n, data1, data2) / norm

  mean2 <- triangular.int(
    function(theta1,theta2,data1,data2) 
      theta2 * likelihood(theta1,theta2,data1,data2),
    low, high, n, data1, data2) / norm
  
  c(mean1,mean2)
}


# TRAPEZOIDAL INTEGRATION FUNCTION.  Integrates f(x,...) from x=a to x=b,
# using the trapezoidal rule with n+1 points.

trapezoidal.rule <- function(f,a,b,n,...)
{ 
  h <- (b-a) / n
  s <- 0

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

  s*h + (f(a,...)+f(b,...))*(h/2)
}


# TRIANGULAR INTEGRATION FUNCTION.  Integrates f(x,y,...) for y from a to b
# and x from y to b, using n+1 points in each application of the trapeziod rule.

triangular.int <- function (f,a,b,n,...)
{
  trapezoidal.rule (
   function (y,f,b,n,...) trapezoidal.rule(f,y,b,n,y,...),
   a, b, n, f, b, n, ...)
}


# DATA TO TEST ON.

data1 <- c(6.1, 7.2, 6.9, 5.4)
data2 <- c(6.2, 5.7, 6.4, 5.2)