# Homework 4

library('MCMCpack')
N=20
y=matrix(ncol=2,nrow=N)
y[,1]=rnorm(N,10,1)
y[,2]=rnorm(N,10.5,1)

#set prior parameters
theta1=0
theta2=0
delta1=1000
delta2=1000
a=.1
b=.1

#stuff for chain
M=1000
mu1=1:M
mu2=1:M
s2=1:M
mu1[1]=0
mu2[1]=0
s2[1]=.1

m1=mean(y[,1])
m2=mean(y[,2])

for (m in 2:M)
{
var.prime = 1/((N/s2[m - 1]) + (1/delta1))
mean.prime = ((N * m1/s2[m - 1]) + (theta1/delta1)) * var.prime
mu1[m] = rnorm(1, mean.prime, sqrt(var.prime))
var.prime = 1/((N/s2[m - 1]) + (1/delta2))
mean.prime = ((N * m2/s2[m - 1]) + (theta2/delta2)) * var.prime
mu2[m] = rnorm(1, mean.prime, sqrt(var.prime))
s2.scale = sum((y[, 1] - mu1[m])^2)/2 + sum((y[, 2] - mu2[m])^2)/2 + b
s2[m] = rinvgamma(1, shape = a + N, scale = s2.scale)
}

hist(mu1-mu2)

#################################
#New Analysis For Alpha Model:
#Consider the model: #y1~Normal(mu+alpha,sigma^2) #y2~Normal(mu-alpha,sigma^2)

theta=0
delta=1

#stuff for chain
M=1000
alpha=1:M
mu=1:M
s2=1:M
alpha[1]=0
mu[1]=10
s2[1]=.1

for (m in 2:M)
{
#alpha
var.prime = 1/((2*N/s2[m - 1]) + (1/delta))
mean.prime = (N*(m1-m2)/s2[m-1]) * var.prime
alpha[m] = rnorm(1, mean.prime, sqrt(var.prime))
#mu
var.prime = s2[m - 1]/(2*N)
mean.prime = (m1+m2)/2
mu[m] = rnorm(1, mean.prime, sqrt(var.prime))
#s2
s2.scale = sum((y[,1]-mu[m]-alpha[m])^2)/2 + sum((y[,2]-mu[m]+alpha[m])^2)/2 + b
s2[m] = rinvgamma(1, shape = a + N, scale = s2.scale)
}

plot(density(mu1-mu2),col='darkgreen',ylim=c(0,1.5))
lines(density(2*alpha),col='darkred')