#######################################
#	samples mu, conditional on sdev
#######################################
samplemu = function(xbar,sdev,N){
	return(rnorm(1,xbar,sdev/sqrt(N)))
}

#######################################
# samples sdev, conditional on mu
# This is actually a draw from an
# inverse chi distribution
#######################################
samplesd = function(X,mu,N){
	sqrt(sum((X-mu)^2)/rchisq(1,N))
}

#######################################
#	samples mu and sdev; returns list
#	consisting of the samples
#######################################
samplen = function(n=1000,mu=0, sdev=1) {
	xbar = mean(X)
	N = length(X)
	
 	mus = rep(NaN,n) #Create vectors to hold data
	sdevs = rep(NaN,n)
	
	for(i in 1:n) {
		mus[i] = mu = samplemu(xbar,sdev,N)
		sdevs[i] = sdev = samplesd(X,mu,N)
	}
	list(mu=mus,sd=sdevs) #return results in a list
}

#######################################
#	generate data
#######################################
X = rnorm(15,3,5) #mean 3, sd 5

#######################################
#	run program and display results
#######################################
run = function(n=1000,mu=0, sdev=1) {
	z=samplen(n,mu,sdev)
	quartz("muhist")
	hist(z$mu) #mu marginal
	quartz("sdhist")
	hist(z$sd) #sd marginal
	quartz("muplot")
	plot(z$mu) #look at sequence
	quartz("sdplot")
	plot(z$sd) #look at sequence
	cat("mubar = ",mean(z$mu),"\n") 	#posterior mean of mu
	cat("sdbar = ",mean(z$sd),"\n") 	#posterior mean of sd
	quartz("correlation plot") #look at 	correlations
	plot(unclass(z$mu),unclass(z$sd))
}

run()