############################################
#	This function plots a histogram for 
#  objects in class samp
#
############################################

hist.samp = function(vals, varname="Variable", adj=2)
{
	x=unclass(vals)
	quartz(varname)
	hist(x,50,freq=F,main=paste("Histogram of", varname), xlab=varname)
	lines(density(x,adjust=adj))
	rug(x)
}

#############################################	
#	This function plots the time history 
#  of a sample in class samp
#
############################################

plot.samp = function(vals,varname="x")
{
	x=unclass(vals)
	quartz(varname)
	plot(x,xlab="Sample Number",ylab=varname)
}

#############################################	
#	This function plots the time history 
#  and histogram of the first and last
#	elements of the array of samples
#	of theta
#
############################################
plotthetas = function(x){
	N = dim(x)[2]
	quartz("thetaplot",width=8,height=8)
	par(mfrow=c(2,2))
	plot(x[,1])
	hist(x[,1],50)
	plot(x[,N])
	hist(x[,N],50)
}

# Generate data
X=c(0.4,0.378,0.356,0.333,0.311,0.311,0.289,
  0.267,0.244,0.244,0.222,0.222,0.222,0.222,
  0.222,0.200,0.178,0.156)
	sigma = sqrt(0.25*(1-0.25)/45)
N = length(X)

samplen = function(n=1000){
	theta = X #Starting values
	mu = mean(X)
	tau = sd(X)
	mus = rep(NaN,n)
	taus = rep(NaN,n)
	thetas = matrix(NaN,n,N)
	
	for(i in 1:n){
		taus[i] = tau = 
			sqrt(sum((mu-theta)^2)/rchisq(1,N-1))
		mus[i] = mu = 
			rnorm(1,mean(theta),tau/sqrt(N))
		S = 1/(1/sigma^2+1/tau^2) 
		a = (X/sigma^2+mu/tau^2)*S
		thetas[i,1:N] = theta = rnorm(N,a,sqrt(S))
	}
	class(taus) = "samp"
	class(mus) = "samp"
	list(tau=taus,mu=mus,theta=thetas)
}

#######################################
#	run program and display results
#######################################
run = function(n=1000){
	z = samplen(n)
	plot(z$mu)
	hist(z$mu)
	plot(z$tau)
	hist(z$tau)
	plotthetas(z$theta)
	cat("mubar  = ",mean(z$mu),"\n") 	#posterior mean of mu
	cat("taubar = ",mean(z$tau),"\n") 	#posterior mean of sd
	cat("theta1bar = ",mean(z$theta[,1]),"\n")
	cat("thetaNbar = ",mean(z$theta[,N]),"\n")
}

run()