#################
# PART I - ii)
#################
a     = c(2,2,3,2)
b     = c(1,2,2,3)
theta = seq(0,10,length=1000)

par(mfrow=c(1,1))
plot(theta,dgamma(theta,a[1],b[1]),xlab=expression(theta),
       type="l",ylab="Probability density function (pdf)",
       main="",lwd=2,ylim=c(0,1.2))
for (i in 2:4)
  lines(theta,dgamma(theta,a[i],b[i]),col=i,lwd=2)
legend("topright",legend=c(
paste("Gamma(",a[1],",",b[1],")",sep=""),
paste("Gamma(",a[2],",",b[2],")",sep=""),
paste("Gamma(",a[3],",",b[3],")",sep=""),
paste("Gamma(",a[4],",",b[4],")",sep="")),col=1:4,lwd=2,bty="n")

mu = a/b
sig2 = a/b^2
cbind(a,b,mu,sig2,sqrt(sig2),qgamma(0.025,a,b),qgamma(0.975,a,b))
# a b        mu      sig2       sig        95% CI
# 2 1 2.0000000 2.0000000 1.4142136 0.24220928 5.571643
# 2 2 1.0000000 0.5000000 0.7071068 0.12110464 2.785822
# 3 2 1.5000000 0.7500000 0.8660254 0.30933606 3.612344
# 2 3 0.6666667 0.2222222 0.4714045 0.08073643 1.857214

#################
# PART I - iii)
#################
Ms    = 10^(3:5)
M     = 10^5
draws = matrix(0,M,4)
set.seed(12345)
for (j in 1:4)
  draws[,j] = rgamma(M,a[j],b[j])

par(mfrow=c(3,4))
for (i in 1:3)
  for (j in 1:4){
    hist(draws[1:Ms[i],j],prob=TRUE,main="",xlab=paste("M=",Ms[i],sep=""))
    lines(theta,dgamma(theta,a[j],b[j]),col=2,lwd=2)
    legend("topright",legend=paste("Gamma(",a[j],",",b[j],")",sep=""),
           bty="n")
}
mus   = matrix(0,3,4)
sig2s = matrix(0,3,4)
for (i in 1:3)
for (j in 1:4){
  mus[i,j]   = mean(draws[1:Ms[i],j])
  sig2s[i,j] = var(draws[1:Ms[i],j])
}
mus   = rbind(mus,mu)
sig2s = rbind(sig2s,sig2)
rownames(mus) = c("1000","10000","1000000","Exact")
rownames(sig2s) = c("1000","10000","1000000","Exact")
round(mus,4)
round(sig2s,4)

#################
# PART I - iv)
#################
R     = 100
M     = 10^5
draws = array(0,c(R,M,4))
set.seed(12345)
for (r in 1:R)
  for (j in 1:4)
    draws[,,j] = rgamma(R*M,a[j],b[j])

mus1   = array(0,c(R,3,4))
sig2s1 = array(0,c(R,3,4))
for (i in 1:R)
  for (j in 1:3){
    mus1[i,j,]   = apply(draws[i,1:Ms[j],],2,mean)
    sig2s1[i,j,] = apply(draws[i,1:Ms[j],],2,var)
  }
  
par(mfrow=c(2,4))
for (i in 1:4){
  boxplot(mus1[,,i],names=Ms,xlab="MC size",ylab="Gamma mean")
  abline(h=mu[i],col=2)
  title(paste("Gamma(",a[i],",",b[i],")",sep=""))
}
for (i in 1:4){
  boxplot(sig2s1[,,i],names=Ms,xlab="MC size",ylab="Gamma variance")
  abline(h=sig2[i],col=2)
}

#################
# PART II - i)
#################
nu = 3
M  = 5000
N  = 1000

set.seed(31415926)
theta = rt(M,df=nu)
w.a = exp(-0.5*theta^2)/(1+theta^2/nu)^(-(nu+1)/2)
w.b = exp(-abs(theta))/(1+theta^2/nu)^(-(nu+1)/2)
theta.a = sample(theta,size=N,replace=TRUE,prob=w.a)
theta.b = sample(theta,size=N,replace=TRUE,prob=w.b)

th.range = range(c(theta,theta.a,theta.b))
thetas = seq(th.range[1],th.range[2],length=1000)

par(mfrow=c(1,3))
hist(theta,prob=TRUE,main="Proposal draws",xlab=expression(theta))
lines(thetas,dt(thetas,df=nu),col=2)
hist(theta.a,prob=TRUE,main="Draws from target A",xlab=expression(theta))
lines(thetas,(2*pi)^(-0.5)*exp(-0.5*thetas^2),col=2)
hist(theta.b,prob=TRUE,main="Draws from target B",xlab=expression(theta))
lines(thetas,exp(-abs(thetas))/2,col=2)

#################
# PART II - ii)
#################
c(mean(theta.a),var(theta.a))
c(mean(theta.b),var(theta.b))


#################
# PART II - iii)
#################
M  = 50000
N  = 10000

set.seed(3141592)
theta = runif(M,-10,10)
w.a = exp(-0.5*theta^2)
w.b = exp(-abs(theta))
theta.a = sample(theta,size=N,replace=TRUE,prob=w.a)
theta.b = sample(theta,size=N,replace=TRUE,prob=w.b)

th.range = range(c(theta,theta.a,theta.b))
thetas = seq(th.range[1],th.range[2],length=1000)

par(mfrow=c(1,3))
hist(theta,prob=TRUE,main="Proposal draws",xlab=expression(theta))
abline(h=1/20,col=2)
hist(theta.a,prob=TRUE,main="Draws from target A",xlab=expression(theta))
lines(thetas,(2*pi)^(-0.5)*exp(-0.5*thetas^2),col=2)
hist(theta.b,prob=TRUE,main="Draws from target B",xlab=expression(theta))
lines(thetas,exp(-abs(thetas))/2,col=2)

c(mean(theta.a),var(theta.a))
c(mean(theta.b),var(theta.b))



#################
# PART II - iv)
#################
tail.a = 1-pnorm(2)
x      = seq(0,2,length=100000)
h      = x[2]-x[1]
tail.b = 0.5-sum(0.5*h*exp(-x))
c(tail.a,tail.b)

M  = 100000
N  = 20000


nu = 3
set.seed(31415926)
theta = rt(M,df=nu)
w.a = exp(-0.5*theta^2)/(1+theta^2/nu)^(-(nu+1)/2)
w.b = exp(-abs(theta))/(1+theta^2/nu)^(-(nu+1)/2)
theta.a1 = sample(theta,size=N,replace=TRUE,prob=w.a)
theta.b1 = sample(theta,size=N,replace=TRUE,prob=w.b)


set.seed(314159)
theta = runif(M,-10,10)
w.a = exp(-0.5*theta^2)
w.b = exp(-abs(theta))
theta.a2 = sample(theta,size=N,replace=TRUE,prob=w.a)
theta.b2 = sample(theta,size=N,replace=TRUE,prob=w.b)

c(tail.a,mean(theta.a1>2),mean(theta.a2>2))
c(tail.b,mean(theta.b1>2),mean(theta.b2>2))




# Replication
R  = 1000
M  = 10000
N  = 2000
nu = 3
set.seed(3145)
tails = matrix(0,R,4)
for (r in 1:R){
  theta = rt(M,df=nu)
  w.a = exp(-0.5*theta^2)/(1+theta^2/nu)^(-(nu+1)/2)
  w.b = exp(-abs(theta))/(1+theta^2/nu)^(-(nu+1)/2)
  theta.a1 = sample(theta,size=N,replace=TRUE,prob=w.a)
  theta.b1 = sample(theta,size=N,replace=TRUE,prob=w.b)

  theta = runif(M,-10,10)
  w.a = exp(-0.5*theta^2)
  w.b = exp(-abs(theta))
  theta.a2 = sample(theta,size=N,replace=TRUE,prob=w.a)
  theta.b2 = sample(theta,size=N,replace=TRUE,prob=w.b)

  tails[r,1:2] = c(mean(theta.a1>2),mean(theta.a2>2))
  tails[r,3:4] = c(mean(theta.b1>2),mean(theta.b2>2))
}

par(mfrow=c(1,2))
boxplot(tails[,1:2])
boxplot(tails[,3:4])