3.4 Bayesian Robustness (pp. 44-49).txt

#An Illustration of Bayesian Robustness
#(pp.44-49)
#prior belief - Joe's IQ median is 100 and 90% sure it is between 80 to 120
#we use normal.select function in LearnBayes to find 
#the values of the mean and standard deviation of the normal density that match
#our prior beliefs

library(LearnBayes)
quantile1=list(p=.5,x=100); quantile2=list(p=.95,x=120)
normal.select(quantile1,quantile2)

#prior is normal with median=100 and 95%=120
#see derivation of posterioir mean and sigma in p.46

mu=100
tau=12.16
sigma=15
n=4
se=sigma/sqrt(4)
ybar=c(110,125,140)
tau1=1/sqrt(1/se^2+1/tau^2)
mu1=(ybar/se^2+mu/tau^2)*tau1^2
summ1=cbind(ybar, mu1, tau1)
summ1

#we try a new prior, t density (median=100, scale t=? and 2 df)
#get scale t. See derivation in p.47
#we find a scale t such that the 95% is 120

tscale=20/qt(.95,2)
tscale

#display the normal and t priors in a single graph
#dt returns the height of the probability density function

par(mfrow=c(1,1))
curve(1/tscale*dt((x-mu)/tscale,2), 
from=60, to=140, xlab="theta", ylab="Prior Density")
curve(dnorm(x,mean=mu,sd=tau), add=TRUE, lwd=3)
legend("topright",legend=c("t density","normal density"), 
lwd=c(1,3))

#We use "prior time likelihood approach"
#Here is the algorithm
#http://www.cyclismo.org/tutorial/R/probability.html
#http://seankross.com/notes/dpqr/
#dnorm give the height of any x (i.e. pdf of normal)

norm.t.compute=function(ybar){
theta=seq(60,180,lenth=500)
like=dnorm(theta,mean=ybar,sd=sigma/sqrt(n))
prior=dt((theta-mu)/tscale,2)
post=prior*like
post=post/sum(post)
m=sum(theta*post)
s=sqrt(sum(theta^2*post)-m^2)
c(ybar, m, s)}

summ2=t(sapply(c(110, 125, 140), norm.t.compute))
dimnames(summ2)[[2]]=c("ybar", "mu1 t", "tau1 t")
summ2

#combine posterior moments of theta using two priors by 
#combining the two R matrices summ1 and summ2

cbind(summ1,summ2)

#Graphs the posterior densities for two priors

theta=seq(60, 180, length=500)
normpost=dnorm(theta, mu1[3], tau1)
normpost=normpost/sum(normpost)
plot(theta, normpost, type="l", lwd=3,ylab="Posterior Density")
like=dnorm(theta,mean=140,sd=sigma/sqrt(n))
prior=dt((theta-mu)/tscale,2)
tpost=prior*like/sum(prior*like)
lines(theta,tpost)
legend("topright", legend=c("t prior", "normal prior"), lwd=c(1,3))