make.prior <- function(a,b,log=FALSE) {
  function(theta) {
    if ( length(theta) > 1 ) stop("theta must by a scalar.")
    dgamma(theta,a,b,log=log)
  }
}

make.likelihood <- function(data,log=FALSE) {
  function(theta) {
    if ( length(theta) > 1 ) stop("theta must by a scalar.")
    f <- function(x) {
      dnorm(x,mean=theta,sd=sqrt(theta),log=log)
    }
    ifelse(log,sum(sapply(data,f)),prod(sapply(data,f)))
  }
}

make.posterior <- function(prior,likelihood,normalizer=ifelse(get("log",environment(prior)),0,1),log=FALSE) {
  prior.on.log <- get("log",environment(prior))
  if ( get("log",environment(likelihood)) != prior.on.log ) stop("The scales of the two functions are not the same.")
  if ( log ) {
    if ( prior.on.log ) function(theta) { likelihood(theta) + prior(theta) + normalizer }
    else function(theta) { log(likelihood(theta)) + log(prior(theta)) + log(normalizer) }
  } else {
    if ( prior.on.log ) function(theta) { exp(likelihood(theta) + prior(theta) + normalizer) }
    else function(theta) { likelihood(theta) * prior(theta) * normalizer }
  }
}

integrate <- function(f,method="trapezoidal.rule",xlim,M) {
  if ( method == "trapezoidal.rule" ) {
    h <- (xlim[2]-xlim[1])/M
    x <- seq(xlim[1],xlim[2],by=h)
    y <- sapply(x,f)
    h*((y[1]+y[length(y)])/2 + sum(y[2:(length(y)-1)]))
  } else if ( method == "simpsons.rule" ) {
    h <- (xlim[2]-xlim[1])/(2*M)
    x <- seq(xlim[1],xlim[2],by=h)
    y <- sapply(x,f)
    h/3 * sum( c(1,rep(c(4,2),times=M-1),4,1) * y )
  } else if ( method == "uniform" ) {
    x <- runif(M,xlim[1],xlim[2])
    y <- sapply(x,f) / ( 1 / (xlim[2]-xlim[1]) )
    mean(y)
  } else stop("Method not recognized")
}



data <-  c(3.08, 0.68, 2.09, 0.87, -0.02, 0.25, 1.98, 1.47, 1.95, 0.99)

x <- seq(0,1.2*max(data),length=1000)

prior <- make.prior(2,2,log=TRUE)
prior <- make.prior(0.001,0.001,log=TRUE)
prior <- make.prior(1,1,log=TRUE)
y.prior <- exp(sapply(x,prior))
plot(x,y.prior,type="l",xlab=expression(theta),ylab=expression(paste("p(",theta,")")),main="Prior Density")

likelihood <- make.likelihood(data,log=TRUE)
y.likelihood <- exp(sapply(x,likelihood))
plot(x,y.likelihood,type="l",xlab=expression(theta),ylab=expression(paste("p(",theta,")")),main="Likelihood Function")

unnormalized.posterior <- make.posterior(prior,likelihood,log=FALSE)
y.posterior <- sapply(x,unnormalized.posterior)
plot(x,y.posterior,type="l",xlab=expression(theta),ylab=expression(paste("p(",theta,")")),main="Unnormalized Posterior Density")

est <- (mean(data)+var(data))/2
range <- est+c(-2,3)*sqrt(est/length(data))
range

system.time(normalizer <- integrate(unnormalized.posterior,"uniform",range,1000)); normalizer
system.time(normalizer <- integrate(unnormalized.posterior,"trapezoidal.rule",range,100)); normalizer
system.time(normalizer <- integrate(unnormalized.posterior,"simpsons.rule",range,50)); normalizer
posterior <- make.posterior(prior,likelihood,-log(normalizer),log=FALSE)
y.posterior <- sapply(x,posterior)
plot(x,y.posterior,type="l",xlab=expression(theta),ylab=expression(paste("p(",theta,")")),main="Prior and Posterior Density")
lines(x,y.prior,col="red")

one <- integrate(posterior,"trapezoidal.rule",range,1000)
one

mean <- integrate(function(theta) theta * posterior(theta),"trapezoidal.rule",range,1000)
mean

mode <- nlm(function(theta) -posterior(theta),mean(data))$estimate
mode

mle1 <- nlm(function(theta) -likelihood(theta),mean(data))$estimate
mle1

mle2 <- optimize(function(theta) -likelihood(theta),c(range[1],range[2]))$minimum
mle2