# This is a R script file for the one sample normal problem. 


##############################################################################
# Open a graphics window

X11()
par(mfrow=c(2,2))

##############################################################################
# Generate the simulated data.  Take a random sample of size n from a
# Normal Distribution with mean mu and standard deviations sigma.

n <- 30
mu <- 10
sigma <- 1
sigmasq <- sigma**2
tau <- 1/sigmasq

y <- rnorm(n,mu,sigma)

y.bar <- mean(y)
y.var <- var(y)
y.sd <- sqrt(y.var)

y.min <- floor(min(y))		# This gives the largest integer less than 
				# all of the y's.
y.max <- ceiling(max(y))	# This gives the smallest integer greater
				# than all of the y's.
br <- seq(y.min,y.max,0.5) 
hist(y, desity=-1, breaks = br, main="histogram of the data, y")

##############################################################################
# Command to use the sta145 splus library

library(sta145, lib.loc='/pkg/splus-library')

##############################################################################
# Plot the priors.
# Prior on mu.

# We think that mu is about 8 and we are 95% sure it is between 5 and 11.
# After playing around a bit I came up with the following. 

m0 <- 8
s0 <- 1.53
s0sq <- s0**2
t0 <- 1/s0sq

mup.05 <- qnorm(0.025,m0,s0)
mup.05
mup.95 <- qnorm(0.975,m0,s0)
mup.95

x <- seq((m0-4*s0),(m0+4*s0),0.01)
d.mu <- dnorm(x,m0,s0)
plot(x,d.mu,type="l",main="prior on mu")

# Prior on tau.

# We think sigma is about 1 and we are 95% sure it is between 0.25 and
# 1.75.  This implies sigmasq is between (0.25)^2 and (1.75)^2.  And
# finally, tau is between 1/(1.75)^2 =  0.33 and 1/(0.25)^2 = 16.
# After playing around a bit the closest I could get was the following.

# Note the mode of a gamma is a/(b-1)

mode <- 1

a <- 5.5
b <- (a - 1)/mode

g1 <- a/2
g2 <- b/2

taup.05 <- qgamma(0.05,g1,g2)
taup.05
taup.95 <- qgamma(0.95,g1,g2)
taup.95

x <- seq(0,( g1/g2 + 4*sqrt(g1/g2^2) ), 0.001)
d.tau <- dgamma(x,g1,g2)
plot(x,d.tau,type="l",main="prior on tau")

##############################################################################
# Composition Rule: Generate a random sample of ( mu(i),tau(i) )

M <- 1000
BurnIn <- 100 

mu.sample <- numeric(M)
tau.sample <- numeric(M)

start <- 1
tau.in <- start

for (i in 1:M) {
	mu.star <-  (n*tau.in/(n*tau.in+t0))*y.bar + (t0/(n*tau.in+t0))*m0
	tau.star <- n*tau.in + t0
	mu.sample[i] <- rnorm(1,mu.star,sqrt(1/tau.star))
	mu.in <- mu.sample[i]
	g1.star <- (n+a)/2
	g2.star <- (n*(y.bar-mu.in)**2 + (n-1)*y.var + b)/2
	tau.sample[i] <- rgamma(1,g1.star,g2.star)
	tau.in <- tau.sample[i]
	#cat(i)
}

##############################################################################
# Posteriors

X11()
par(mfrow=c(2,2))

# Crude estimate of p(mu|Y)

mu.min <- floor(min(mu.sample[BurnIn:M]))
mu.max <- ceiling(max(mu.sample[BurnIn:M]))
br <- seq(mu.min,mu.max,0.05)   
hist(mu.sample[BurnIn:M], probability = T, density=-1, breaks = br,
main="posterior of mu")

# Crude estimate of p(tau|Y).

tau.min <- floor(min(tau.sample[BurnIn:M]))
tau.max <- ceiling(max(tau.sample[BurnIn:M]))
br <- seq(tau.min,tau.max,0.05)
hist(tau.sample[BurnIn:M], probability = T, density=-1, breaks = br,
main="posterior of tau")

# Crude estimate of mu.

mu.pm <- mean(mu.sample[BurnIn:M])
mu.pv <- var(mu.sample[BurnIn:M])
mu.psd <- sqrt(var(mu.sample[BurnIn:M]))
mu.pCI <- quantile(mu.sample[BurnIn:M],c(0.05,0.95))

# Crude estimate of tau.

tau.pm <- mean(tau.sample[BurnIn:M])
tau.pv <- var(tau.sample[BurnIn:M])
tau.psd <- sqrt(var(tau.sample[BurnIn:M]))
tau.pCI <- quantile(tau.sample[BurnIn:M],c(0.05,0.95)) 

# Crude estimate of sigmasq.

sigmasq.pm <- mean(1/tau.sample[BurnIn:M])
sigmasq.pv <- var(1/tau.sample[BurnIn:M])
sigmasq.psd <- sqrt(var(1/tau.sample[BurnIn:M]))
sigmasq.pCI <- quantile((1/tau.sample[BurnIn:M]),c(0.05,0.95)) 

# Crude estimate of sigma.

sigma.pm <- mean(1/sqrt(tau.sample[BurnIn:M]))
sigma.pv <- var(1/sqrt(tau.sample[BurnIn:M]))
sigma.psd <-sqrt(var(1/sqrt(tau.sample[BurnIn:M])))
sigma.pCI <- quantile((1/sqrt(tau.sample[BurnIn:M])),c(0.05,0.95)) 

# Better posterior estimate of mu, using muhat(tau).
# Better posterior estimate of p(mu|Y).

# Note: mu.star here is muhat(tau) in class.

mu.star <- ( n*tau.sample[(BurnIn+1):M]/
		(n*tau.sample[(BurnIn+1):M] + t0))*y.bar
                + ( t0/(n*tau.sample[(BurnIn+1):M] + t0) )*m0
muhat <- mu.star
mu.better <- mean(muhat)

tau.star <- (n*tau.sample[(BurnIn+1):M] + t0)

x <- seq(8.5,11,0.01)
px <- matrix(0,(M-BurnIn),length(x))
for (i in 1:(M-BurnIn))	{
	px[i,] <- dnorm(x,mu.star[i],sqrt((1/tau.star[i])))
}
px <- apply(px,2,mean)
X11()
par(mfrow=c(2,2))
plot(x,px,type='l', main="Better estimate of p(mu|Y)")

# Better posterior estimate of tau, using tauhat(mu).
# Better posterior estimate of p(tau|Y). 

g1.star <- (n+a)/2
g2.star <- (n*((y.bar-mu.sample[(BurnIn+1):M])**2)+(n-1)*y.var+b)/2
tauhat <- g1.star/g2.star
tau.better <- mean(tauhat)

x <- seq(0,2,0.01)
px <- matrix(0,(M-BurnIn),length(x))
for (i in 1:(M-BurnIn)) {
        px[i,] <- dgamma(x,g1.star,g2.star[i]) 
}
px <- apply(px,2,mean)
plot(x,px,type='l',main="Better estimate of p(tau|Y)")

#####################################################~#########################
# The Markov Chains

X11()
par(mfrow=c(2,2))
ts.plot(mu.sample[BurnIn:M], main="Gibbs MC Chain for mu, BurnIn:M")
acf(mu.sample[BurnIn:M])
ts.plot(tau.sample[BurnIn:M],main="Gibbs MC Chain for tau, BurnIn:M")
acf(tau.sample[BurnIn:M])

# Monitor convergence

muMean <- numeric(M)
muSd <- numeric(M)
tauMean <- numeric(M)
tauSd <- numeric(M)
for (i in 1:M) {
	muMean[i] <- mean(mu.sample[1:i])
	muSd[i] <- sqrt(var(mu.sample[1:i]))
	tauMean[i] <- mean(tau.sample[1:i])
	tauSd[i] <- sqrt(var(tau.sample[1:i]))
}

X11()
par(mfrow=c(2,2))
ts.plot(muMean, main="Cumulative Mean of Gibbs MC Chain for mu")
ts.plot(muSd, main="Cumulative Sd of Gibbs MC Chain for mu")      
ts.plot(tauMean, main="Cumulative Mean of Gibbs MC Chain for tau") 
ts.plot(tauSd, main="Cumulative Sd of Gibbs MC Chain for tau")

##############################################################################
# Prediction

Ynew <- rnorm((M-BurnIn),mu.sample[BurnIn:M],sqrt(tau.sample[BurnIn:M]))

Ynew.min <- floor(min(Ynew))
Ynew.max <- ceiling(max(Ynew))
br <- seq(Ynew.min,Ynew.max,0.5)
X11()
hist(Ynew, probability=T, density=-1,breaks=br,main="predicitve
distribution")

# Crude extimate of Ynew.

Ynew.m <- mean(Ynew)
Ynew.v <- var(Ynew)
Ynew.sd <- sqrt(var(Ynew))
Ynew.CI <- quantile(Ynew, c(0.05,0.95))

##############################################################################