###############################################
 #Prior Sensitivity Example
 #Temperature Data
 #1/12/06
 ###############################################

#Set working directory
setwd("U:\\Teaching\\Stat694\\Computing\\")

#Read in the temperature data
temp.data <- read.table("temperature.dat",header=T)
temp <- temp.data$Temp
N <- length(temp)

#Set up model
sigma.y <- 2          #assumed standard deviation of the data
prior.mean <- 26      #prior mean
prior.sd <- 2         #prior standard deviation

#Examine how the posterior distribution changes as the number of data points increases
par(mfrow=c(1,1),ask=T)
for(n in 1:N){
      y <- temp[1:n]    #select the first n measurements

      post.sd <- sqrt(1/((1/prior.sd^2)+(n/sigma.y^2)))     #calculate posterior sd
      post.mean <- post.sd^2*((prior.mean/prior.sd^2) + (n*mean(y)/sigma.y^2))  #calculate posterior mean

      #construct plots
      x <- seq(15,50,length=500)
      plot(x,dnorm(x,prior.mean,prior.sd),type='l',col="red",ylim=c(0,dnorm(post.mean,post.mean,post.sd)),xlim = range(x.temp),xlab="temperature",ylab="",main=paste("n = ",n))
      lines(x.temp,dnorm(x.temp,post.mean,post.sd),col="blue")
      points(y,rep(0,n),col="green",pch=19)
      abline(v=prior.mean,col="red")
      abline(v=post.mean,col="blue")
      abline(v=mean(y),col="green")
      legend(40,dnorm(post.mean,post.mean,post.sd),c("prior","data","posterior"),col=c("red","green","blue"),text.col="black",lty=1)
}