# Modification of MDP code to study mixing at a more basic level # First, a Gibbs sampler without the remixing step. If this sampler # is run with a large value of M (so that there are nearly as many # clusters as observations in the data set, mixing is nearly as good # as with the sampler with the remixing step. If the value of M # is relatively small, there are lots of "repeats", as clusters don't # often move. This behavior is not uniform across the cases in the # data set, but depends on the local properties of data and prior. # It is often very different for observations in the middle of the # data set and those at the extremes. # The functions below make use of the quicker code in functions that # end with ".2". These are available in the modified ################################################################ # # Function 5. One iterate of the Gibbs sampler # # Tasks: Generate each theta_i in turn # Generate theta_star # Generate mu # ################################################################ fn.one.iterate.2.nomix <- function(clust,prior,mu,data) { for (i in 1:data$n) { clust <- fn.remove.theta.i(i,clust) clust <- fn.gen.theta.i.2(i,clust,prior,mu,data) } # clust <- fn.gen.theta.star.2(clust,prior,mu,data) mu <- fn.gen.mu.2(clust,prior) ret.obj <- NULL ret.obj$clust <- clust ret.obj$prior <- prior ret.obj$mu <- mu return(ret.obj) } ################################################################ # # Function 6. A brief Gibbs sampler # # Tasks: Set up object (here, matrix) to store results # Run one iterate of Gibbs sampler # Tally results # # Improvements for you to make: # Burn-in -- allow explicit description of burn-in to be # discarded # Subsampling -- Not to be done unless storage issues are # important. But, allow subsampling of the # output # Initialization -- an automated initialization routine. # Best to allow a couple of options for # the initialization. # ################################################################ fn.gibbs.sampler.2.nomix <- function(n.reps,prior,data,clust,mu) { # Insert initialization routine if desired # Insert burn-in period if desired res.mat <- matrix(rep(0,n.reps*(data$n+1)),nrow=n.reps) for (i in 1:n.reps) { tmp <- fn.one.iterate.2.nomix(clust,prior,mu,data) clust <- tmp$clust mu <- tmp$mu res.mat[i,] <- c(mu,clust$theta.star[clust$s]) # print(clust$k) } return(res.mat) } # create appropriate data object data1 <- NULL data1$x <- bodyfat.mat[,4] data1$n <- length(data1$x) prior1 <- NULL prior1$mu.0 <- 180 prior1$rho.sq <- 400 prior1$M <- 1 # for this data set, use M=20 for a # large value, M=1 for a small value prior1$tau.sq <- 225 prior1$sig.sq <- 100 clust1 <- NULL clust1$k <- data1$n clust1$s <- 1:data1$n clust1$n.i <- rep(1,data1$n) clust1$theta.star <- data1$x mu1 <- prior1$mu.0 n.reps <- 2000 date() system.time( res.mat.2.nm <- fn.gibbs.sampler.2.nomix(n.reps,prior1,data1,clust1,mu1) ,) date() par(mfrow=c(2,1)) # choice of individual 82 (position 83 since mu is stored in position 1) # is arbitrary length(unique(res.mat.2.nm[,83])) plot(res.mat.2.nm[,83]) plot(density(res.mat.2.nm[,83])) date() system.time( res.mat.2 <- fn.gibbs.sampler.2(n.reps,prior1,data1,clust1,mu1) ,) date() length(unique(res.mat.2[,83])) plot(res.mat.2[,83]) plot(density(res.mat.2[,83])) acf(res.mat.2.nm[,83]) acf(res.mat.2[,83]) # Note--the long, horizontal black blobs in the "no-mixing" plot # indicate poor mixing. Also note the density plots. The # auto-correlation plots tell an interesting tale. # Next, a function that lets us make use of a "denoised" version # of theta_i. This denoised version replaces the generated value # of theta_i with the mean of the appropriate theta.star value (just # before its generation). While these "denoised" theta_i would not # be used to get estimates of the posterior density of theta_i or # its variance, they can be used for "Rao-Blackwellized" estimates # of the posterior mean of theta_i. # of the time plot for the theta_i. ################################################################ # # Function 3. A function to generate theta_star # # Tasks: Loop through i = 1, ..., k # Find cond'l posterior distribution for theta_star[i] # Generate theta_star[i] # ################################################################ fn.gen.theta.star.3 <- function(clust,prior,mu,data) { k <- clust$k n.i <- clust$n.i s <- clust$s theta.star <- clust$theta.star tau.sq <- prior$tau.sq sig.sq <- prior$sig.sq x <- data$x ############################# # for (i in 1:k) # { # tmp.m <- ((n.i[i]/sig.sq)*mean(x[s==i]) + # (1/tau.sq)*mu) / # ((n.i[i]/sig.sq) + (1/tau.sq)) # tmp.v <- 1/((n.i[i]/sig.sq) + (1/tau.sq)) # theta.star[i] <- rnorm(n=1,mean=tmp.m,sd=sqrt(tmp.v)) # } ### replacement for above ### tmp.m <- rep(0,k) for (i in 1:k) { tmp.m[i] <- ((n.i[i]/sig.sq)*mean(x[s==i]) + (1/tau.sq)*mu) / ((n.i[i]/sig.sq) + (1/tau.sq)) } tmp.v <- 1/((n.i/sig.sq) + (1/tau.sq)) theta.star <- rnorm(n=k,mean=tmp.m,sd=sqrt(tmp.v)) ############################# clust$theta.star <- theta.star # The single line below passes the means of the theta.star clust$theta.star.mean <- tmp.m return(clust) } ################################################################ # # Function 5. One iterate of the Gibbs sampler # # Tasks: Generate each theta_i in turn # Generate theta_star # Generate mu # ################################################################ fn.one.iterate.3 <- function(clust,prior,mu,data) { for (i in 1:data$n) { clust <- fn.remove.theta.i(i,clust) clust <- fn.gen.theta.i.2(i,clust,prior,mu,data) } clust <- fn.gen.theta.star.3(clust,prior,mu,data) mu <- fn.gen.mu.2(clust,prior) ret.obj <- NULL ret.obj$clust <- clust ret.obj$prior <- prior ret.obj$mu <- mu return(ret.obj) } ################################################################ # # Function 6. A brief Gibbs sampler # # Tasks: Set up object (here, matrix) to store results # Run one iterate of Gibbs sampler # Tally results # # Improvements for you to make: # Burn-in -- allow explicit description of burn-in to be # discarded # Subsampling -- Not to be done unless storage issues are # important. But, allow subsampling of the # output # Initialization -- an automated initialization routine. # Best to allow a couple of options for # the initialization. # ################################################################ fn.gibbs.sampler.3 <- function(n.reps,prior,data,clust,mu) { # Insert initialization routine if desired # Insert burn-in period if desired res.mat <- matrix(rep(0,n.reps*(2*data$n+1)),nrow=n.reps) for (i in 1:n.reps) { tmp <- fn.one.iterate.3(clust,prior,mu,data) clust <- tmp$clust mu <- tmp$mu res.mat[i,] <- c(mu,clust$theta.star[clust$s],clust$theta.star.mean[clust$s]) # print(clust$k) } return(res.mat) } # create appropriate data object data1 <- NULL data1$x <- bodyfat.mat[,4] data1$n <- length(data1$x) prior1 <- NULL prior1$mu.0 <- 180 prior1$rho.sq <- 400 prior1$M <- 20 prior1$tau.sq <- 225 prior1$sig.sq <- 100 clust1 <- NULL clust1$k <- data1$n clust1$s <- 1:data1$n clust1$n.i <- rep(1,data1$n) clust1$theta.star <- data1$x mu1 <- prior1$mu.0 # A run of the code, with 5000 replicates. After the code hsa # been run, a comparison of the variances of the stored values. # Note that the variance of the estimator is not this value # divided by n.reps. An easy way to get a good estimator for # variation in the estimator of theta.i is the batch means method. # See the brief function below. n.reps <- 5000 date() system.time( res.mat.2 <- fn.gibbs.sampler.2(n.reps,prior1,data1,clust1,mu1) ,) date() system.time( res.mat.3 <- fn.gibbs.sampler.3(n.reps,prior1,data1,clust1,mu1) ,) date() var.mat <- matrix(rep(0,252*3),ncol=3) for (i in 1:252) { var.mat[i,1] <- var(res.mat.2[,1+i]) var.mat[i,2] <- var(res.mat.3[,1+i]) var.mat[i,3] <- var(res.mat.3[,253+i]) } apply(var.mat,2,mean) t.test(var.mat[,1] - var.mat[,2]) t.test(var.mat[,1] - var.mat[,3]) t.test(var.mat[,2] - var.mat[,3]) fn.batch.se <- function(x,batchsize=100) { num.mns <- length(x)/batchsize tmp.mns <- rep(0,num.mns) for (i in 1:num.mns) {tmp.mns[i] <- mean(x[((i-1)*batchsize+1):(i*batchsize)])} batch.var <- var(tmp.mns)/num.mns return(sqrt(batch.var)) } mean(apply(res.mat.2[,2:253],2,fn.batch.se)) mean(apply(res.mat.3[,2:253],2,fn.batch.se)) mean(apply(res.mat.3[,254:505],2,fn.batch.se)) par(mfrow=c(2,1)) plot(apply(res.mat.3[,2:253],2,mean),apply(res.mat.3[,2:253],2,var)) plot(apply(res.mat.3[,254:505],2,fn.batch.se),apply(res.mat.3[,2:253],2,fn.batch.se)) lines(apply(res.mat.3[,254:505],2,fn.batch.se),apply(res.mat.3[,254:505],2,fn.batch.se))