R/MixMVN_MCMC.R

In CubicZebra/MVNBayesian: Bayesian Analysis Framework for MVN (Mixture) Distribution

MixMVN_MCMC <-
function(data, steps, pars, values, tol=0.3, random_method = "Fast", ...){
  d <- length(data[1,2][[1]])
  if(any(pars %in% (d+1))){
    stop("Error assignment: conditioned parameters shouldn't include clusters")
  }
  else{
    MCMC <- MixMVN_GibbsSampler(steps, data, random_method = random_method, ...)
    samples_num <- length(MCMC[,1])
    if(length(pars) != length(values) | length(pars) > length(data[1,2][[1]])){
      stop("Error assignment: please confirm that pars and values have the identical and correct dimension  (be no more than that of input data). ")
    }
    else if(length(pars) == 0){
      MCMC <- cbind(MCMC, matrix(data = 0, samples_num, 1))
      MCMC <- cbind(MCMC, matrix(data = 1, samples_num, 1))
    }
    else{
      idx <- matrix(NA, samples_num, 2)
      for (i in 1:samples_num){
        distance <- MCMC[i,pars]-values
        idx[i,1] <- norm(as.matrix(distance), type = "2")
        idx[i,2] <- (idx[i,1] <= tol)
      }
      MCMC <- cbind(MCMC, idx)
    }
    p <- length(MCMC[1,])
    colnames(MCMC)[p-1] <- "Euclidean.Distance"
    colnames(MCMC)[p] <- "Accept"
    AcceptRate <- as.numeric(count(MCMC[,p])[2,2]/samples_num)
    Accept <- MCMC[which(MCMC[,p] == 1),]
    Reject <- MCMC[which(MCMC[,p] == 0),]
    results <- list("AcceptRate"=AcceptRate, "MCMCdata"=MCMC, "Accept"=Accept, "Reject"=Reject)
    return(results)
  }
}