set.seed(11) alpha <- c(0.01, 1.00, 3.00, 5.00) sigma_y <- c(0.50, 1.00, 1.50, 3.00) for (i in 1:length(alpha)) for (j in 1:length(sigma_y)) { results <- crp_gibbs(data = datc, alpha = alpha[i], mu0 = matrix(rep(0, 2), ncol = 2, byrow = TRUE), sigma0 = diag(2) * 3^2, sigma_y = diag(2) * sigma_y[j], c_init = rep(1, nrow(datc))) tab <- apply( results, 1, FUN = function(x) { tab <- table(x) ans <- names(tab[which.max(tab)]) return(ans) }) cat("alpha = ", alpha[i], "sigma_y = ", sigma_y[j], "\n") print(table(tab)) }