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