# run the CRP Gibbs function in Appendix B. alpha <- 0.01 mu0 <- matrix(rep(0, 2), ncol = 2, byrow = TRUE) sigma0 <- diag(2) * 3^2 sigma_y <- diag(2) * 1 c_init <- rep(1, nrow(datc)) results <- crp_gibbs(data = datc, alpha = alpha, mu0 = mu0, sigma0 = sigma0, sigma_y = sigma_y, c_init = rep(1, nrow(datc))) ##### # Gibbs sampling iterations are saved in # `results'. Each row contains the latent class # membership distribution of one person over # MCMC iterations. A person is deemed a member # of c=2 if, e.g., s/he gets assigned most # often to c2. So we tabulate the frequency # counts and whichever cluster with highest # frequency is that person's latent class. ##### tab <- apply( results, 1, FUN = function(x) { tab <- table(x) ans <- names(tab[which.max(tab)]) return(ans) } ) table(tab)