###############################################################
#
# 2nd review session
# Bayesian normal linear regression
#
###############################################################

pdf(file="2nd-reviewsession-graph1.pdf",width=10,height=10)
# National Longitudinal Survey of Youth dataset
data = matrix(scan("logwages-yearseducation.txt"),byrow=TRUE,ncol=2)
n = nrow(data)
y = data[,1]
x = data[,2]

# a)
par(mfrow=c(3,3))
beta0 = 1.17
sig2  = 0.267
b     = 0
B     = 2
beta1 = seq(-7,7,length=10000)
prior = dnorm(beta1,b,sqrt(B))
plot(beta1,prior,type="l",xlab=expression(beta[1]),
     ylab="Density",lwd=3)
title("a) Prior of beta1")

y1   = y-beta0
bhat = sum(x*y1)/(sum(x^2))
beta1 = seq(0.08,0.1,length=10000)
like = dnorm(beta1,bhat,sqrt(sig2/sum(x^2)))
plot(beta1,like,type="l",xlab=expression(beta[1]),
     ylab="Likelihood",lwd=3)
title("a) Likelihood of beta1")

B1 = 1/(1/B+1/sum(x^2))
b1 = B1*(b/B+sum(x*y1)/sum(x^2))
beta1 = seq(0.08,0.1,length=10000)
post = dnorm(beta1,b1,sqrt(B1))
plot(beta1,like,type="l",xlab=expression(beta[1]),
     ylab="Density",lwd=3)
title("a) Posterior of beta1")

# b)
beta1 = 0.091
sig2  = 0.267
b     = 0
B     = 2
beta0 = seq(-7,7,length=10000)
prior = dnorm(beta0,b,sqrt(B))
plot(beta0,prior,type="l",xlab=expression(beta[0]),
     ylab="Density",lwd=3)
title("b) Prior of beta0")

y1   = y-beta1*x
bhat = mean(y1)
beta0 = seq(1.1,1.25,length=10000)
like = dnorm(beta0,bhat,sqrt(sig2/n))
plot(beta0,like,type="l",xlab=expression(beta[0]),
     ylab="Likelihood",lwd=3)
title("b) Likelihood of beta0")

B1 = 1/(n/sig2+1/B)
b1 = B1*(b/B+sum(y1)/sig2)
beta0 = seq(1.1,1.25,length=10000)
post = dnorm(beta0,b1,sqrt(B1))
plot(beta0,like,type="l",xlab=expression(beta[0]),
     ylab="Density",lwd=3)
title("b) Posterior of beta0")

#c)
dig = function(x,a,b){
  dgamma(1/x,a,b)/x^2
}
beta0 = 1.17
beta1 = 0.091
a = 3
b = 0.4
sig2 = seq(0.001,1,length=10000)
prior = dig(sig2,a,b)
plot(sig2,prior,xlab=expression(sig2),type="l",ylab="Density",lwd=3)
title("c) Posterior of sig2")

sig2 = seq(0.2,0.35,length=10000)
like = dig(sig2,(n-2)/2,sum((y-beta0-beta1*x)^2)/2)
plot(sig2,like,xlab=expression(sig2),type="l",ylab="Likelihood",lwd=3)
title("c) Likelihood of sig2")

a1 = a+n/2
b1 = b+sum((y-beta0-beta1*x)^2)/2
sig2 = seq(0.2,0.35,length=10000)
post = dig(sig2,a1,b1)
plot(sig2,post,xlab=expression(sig2),type="l",ylab="Density",lwd=3)
title("c) Posterior of sig2")

dev.off()

pdf(file="2nd-reviewsession-graph2.pdf",width=10,height=10)
#d) 
par(mfrow=c(2,3))
sig2   = 0.267
beta0  = seq(-7,7,length=100)
beta1  = seq(-7,7,length=100)
lprior = matrix(0,100,100)
for (i in 1:100)
  for (j in 1:100)
     lprior[i,j] = dnorm(beta0[i],0,sqrt(2),log=TRUE)+
                   dnorm(beta1[j],0,sqrt(2),log=TRUE)
prior = exp(lprior-max(lprior))

contour(beta0,beta1,prior,xlab=expression(beta[0]),
      ylab=expression(beta[1]),drawlabels=FALSE)
title("d) Prior of (beta0,beta1)")

beta0  = seq(0.9,1.4,length=100)
beta1  = seq(0.06,0.12,length=100)
lprior = matrix(0,100,100)
llike  = matrix(0,100,100)
for (i in 1:100)
  for (j in 1:100)
     llike[i,j] = sum(dnorm(y,beta0[i]+beta1[j]*x,sqrt(sig2),log=TRUE))
     lprior[i,j] = dnorm(beta0[i],0,sqrt(2),log=TRUE)+
                   dnorm(beta1[j],0,sqrt(2),log=TRUE)
prior = exp(lprior-max(lprior))
like = exp(llike-max(llike))

contour(beta0,beta1,like,xlab=expression(beta[0]),
      ylab=expression(beta[1]),drawlabels=FALSE)
title("d) Likelihood of (beta0,beta1)")

lpost = llike + lprior
post  = exp(lpost-max(lpost))

contour(beta0,beta1,post,xlab=expression(beta[0]),
      ylab=expression(beta[1]),drawlabels=FALSE)
title("d) Posterior of (beta0,beta1)")

#e) 
sig2   = 0.267
nu     = 4
beta0  = seq(-7,7,length=100)
beta1  = seq(-7,7,length=100)
lprior = matrix(0,100,100)
for (i in 1:100)
  for (j in 1:100)
     lprior[i,j] = -0.5*(nu+2)*log(1+(1/nu)*(beta0[i]^2+beta1[j]^2))
prior = exp(lprior-max(lprior))

contour(beta0,beta1,prior,xlab=expression(beta[0]),
      ylab=expression(beta[1]),drawlabels=FALSE)
title("e) Prior of (beta0,beta1)")

beta0 = seq(0.9,1.4,length=100)
beta1 = seq(0.06,0.12,length=100)
llike  = matrix(0,100,100)
for (i in 1:100)
  for (j in 1:100){
     llike[i,j] = sum(dnorm(y,beta0[i]+beta1[j]*x,sqrt(sig2),log=TRUE))
     lprior[i,j] = -0.5*(nu+2)*log(1+(1/nu)*(beta0[i]^2+beta1[j]^2))
   }
like  = exp(llike-max(llike))
prior = exp(lprior-max(lprior))

contour(beta0,beta1,like,xlab=expression(beta[0]),
      ylab=expression(beta[1]),drawlabels=FALSE)
title("e) Likelihood of (beta0,beta1)")

lpost = llike + lprior
post  = exp(lpost-max(lpost))

contour(beta0,beta1,post,xlab=expression(beta[0]),
      ylab=expression(beta[1]),drawlabels=FALSE)
title("e) Posterior of (beta0,beta1)")

dev.off()


pdf(file="2nd-reviewsession-graph3.pdf",width=10,height=10)
# f) 
set.seed(1235)
M = 100000

# Proposal/candidate/auxiliary distribution
dmnorm = function(x,m,iS,ldetS){-0.5*ldetS-0.5*t(x-m)%*%iS%*%(x-m)}
X     = cbind(1,x)
bhat  = solve(t(X)%*%X)%*%t(X)%*%y
s2    = mean((y-X%*%bhat)^2)
V     = s2*solve(t(X)%*%X)
iV    = solve(V)
ldetV = det(V,log=TRUE)
betas = matrix(bhat,M,2,byrow=TRUE)+matrix(rnorm(2*M),M,2)%*%chol(V)
w     = rep(0,100000)
for (i in 1:M)
 w[i] = -0.5*(nu+2)*log(1+(1/nu)*(betas[i,1]^2+betas[i,2]^2))+
        sum(dnorm(y,betas[i,1]+betas[i,2]*x,sqrt(sig2),log=TRUE))-
        dmnorm(betas[i,],bhat,iV,ldetV)
w = exp(w-max(w))
w = w/sum(w)
k = sample(1:100000,size=10000,prob=w,replace=TRUE)
betas1 = betas[k,]

par(mfrow=c(1,1))
plot(betas1,xlab=expression(beta[0]),ylab=expression(beta[1]),pch=16)
title("f) Posterior of (beta0,beta1)")
contour(beta0,beta1,post,add=TRUE,drawlabels=FALSE,col=2)

dev.off()

pdf(file="2nd-reviewsession-graph4.pdf",width=10,height=10)
par(mfrow=c(2,2))
# Normal prior
X  = cbind(1,x)
B1 = solve(diag(1/2,2)+t(X)%*%X/sig2)
b1 = B1%*%t(X)%*%y/sig2
E0 = b1[1]+15*b1[2]
V0 = B1[1,1]+225*B1[2,2]+30*B1[1,2]
P  = 1-pnorm(0.1,b1[2],sqrt(B1[2,2]))

# Cauchy prior
E1 = mean(betas1[,2])
V1 = var(betas1[,2])     
P1 = mean(betas1[,2]>0.1)

xxx = seq(min(betas1[,1]),max(betas1[,1]),length=30)
hist(betas1[,1],xlab=expression(beta[0]),ylab="Density",
     prob=TRUE,main="",breaks=xxx)
title("f) Posterior of beta0")
xxx = seq(min(betas1[,1]),max(betas1[,1]),length=1000)
lines(xxx,dnorm(xxx,b1[1],sqrt(B1[1,1])),col=2)
     
xxx = seq(min(betas1[,2]),max(betas1[,2]),length=30)
hist(betas1[,2],xlab=expression(beta[1]),ylab="Density",
     prob=TRUE,main="",breaks=xxx)
title("f) Posterior of beta1")
xxx = seq(min(betas1[,2]),max(betas1[,2]),length=1000)
lines(xxx,dnorm(xxx,b1[2],sqrt(B1[2,2])),col=2)


plot(c(0,10),c(0,10),axes=FALSE,col=0,xlab="",ylab="")
text(5,9.9,"Cauchy prior",cex=1.25)
text(5,9,paste("E(beta1|x,y)     = ",round(E1,3),sep=""),cex=1.25)
text(5,8,paste("V(beta1|x,y)     = ",round(V1,6),sep=""),cex=1.25)
text(5,7,paste("P(beta1>0.1|x,y) = ",round(P1,3),sep=""),cex=1.25)

text(5,5,"Normal prior",cex=1.25,col=2)
text(5,4,paste("E(beta1|x,y)     = ",round(b1[2],3),sep=""),cex=1.25,col=2)
text(5,3,paste("V(beta1|x,y)     = ",round(B1[2,2],6),sep=""),cex=1.25,col=2)
text(5,2,paste("P(beta1>0.1|x,y) = ",round(P,3),sep=""),cex=1.25,col=2)
title("f.1), f.2) e f.3)")

par = betas1[,1]+15*betas1[,2]
xxx = seq(min(par),max(par),length=1000)
hist(par,prob=TRUE,ylab="Density",main="",xlab=expression(beta[0]+15*beta[1]))
title("f.4) Posterior of beta0+15*beta1")
lines(xxx,dnorm(xxx,E0,sqrt(V0)),col=2)

dev.off()