####################################### # 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()