#############################################
# 2.10.
#############################################

  # main functions, need functions defined below.

   # 2.10a
thgrid <- seq(0,1,length=1000)
pgrid  <- post(thgrid)    # note that post works on vector!
                          # alternative use loop "for(i in 1:1000)..."
k <- sum(pgrid*0.001)     # normalization constant
plot(thgrid,pgrid/k,ylim=c(0,1.4),xlab="THETA",ylab="P(TH|Y)", type="l")

   # 2.10b  
th <- rpost(10000)
hist(th,col=0,prob=T)       # compare the histogram with the earlier
                            # plot
   # 2.10c
th <- th[1:1000]            # use the first 1000 only
y6 <- rcauchy(1000,location=th)
hist(y6,col=0,prob=T,xlab="p(Y6|Y1..5)",nbins=20)
                            # note the outliers! The Cauchy has VERY flat tails

  # defines a fucntion which returns the likelihood for th
  # note that f will also work for a whole vector of th's
likl <- function(th)
{
  y <- c(-2,-1,0,1.5,2.5)
  f <- 1
  for (i in 1:5)
    f <- f*1/(1+(y[i]-th)^2)
  f
}

  # defines the prior
prior <- function(th)
{
  1
}

  # posterior
post <- function(th)
{
  prior(th)*likl(th)
}


  # sets up cdf.inv for rpost
  
  # simulating a sample from the posterior
  # note: need to call init.rpost to setup cdf.inv 
  # see book p22/23 for an explaination of the inverse cdf method
rpost <- function(n)
{
  sample(thgrid,size=n,replace=T,prob=pgrid)
}