#Define the posterior assuming y|mu~N(mu,1) & t~Cauchy 
  posterior<-function(mu,y,sigma=1){
     dnorm(mean(y),mu,sigma/sqrt(length(y)))*dt(mu,1)
  }

#Generate a sample with mu=2
  set.seed(2005)
  y<-rnorm(5,2,1)

##################################################################
### Approximate the posterior density, mean and sd on a grid ###
##################################################################

  mu.grid<-seq(0,5,0.01)
  dense<-posterior(mu.grid,y)
  dense.std<-dense/sum(dense)
  mu.mn<-sum(mu.grid*dense.std)
  mu.sd<-sqrt(sum(dense.std*(mu.grid-mu.mn)^2))

  plot(mu.grid,dense.std,type="l",
     main="Approximate posterior distribution",
     xlab=expression(mu),ylab="")
  print("Posterior mean and sd")
  print(round(mu.mn,2))
  print(round(mu.sd,2))

###############################################################################
### Approximate the posterior density, mean and sd using rejection sampling ###
###############################################################################

#pick an enveloping function N(y-bar,s)
  envelope<-function(mu,mn,sd){dnorm(mu,mn,sd)}

#calculate the constant m
  y.bar<-mean(y)
  s<-sd(y)
  dense.env<-envelope(mu.grid,y.bar,s)
  m<-max(dense/dense.env)
  plot(mu.grid,m*dense.env,type="l",lty=2,xlab=expression(mu),ylab="",main="")
  lines(mu.grid,dense)
  legend("topleft",c("Posterior","Envelope"),lty=1:2,inset=.05)

#draw samples
  n.samples<-10000
  mu<-rep(0,n.samples)
  count<-attempts<-0

  while(count<n.samples){
   attempts<-attempts+1
   can<-rnorm(1,y.bar,s)
   prob<-posterior(can,y)/envelope(can,y.bar,s)/m
   if(runif(1)<prob){
      count<-count+1
      mu[count]<-can
   }
  }

  hist(mu,xlab=expression(mu),main="Histogram of the posterior samples")
  print("Posterior mean and sd")
  print(round(mean(mu),2))
  print(round(sd(mu),2))
  print("acceptance ratio")
  print(round(n.samples/attempts,3))