R/sampling_a.R

In BNPMIXcluster: Bayesian Nonparametric Model for Clustering with Mixed Scale Variables

#' @title
#'   MCMC sampling of parameter "\eqn{a}" in the \emph{mixdpcluster} model for Bayesian clustering.
#'
#' @description
#'   Generates a sample from the posterior distribution of \eqn{a} in the \emph{mixdpcluster} model for bayesian clustering.
#'   The simulation is done via Metropolis-Hastings method.
#'
#' @param n number of simulations to generate
#' @param a.ini initialization value
#' @param b parameter \eqn{b} in the posterior distribution of a
#' @param alpha parameter \eqn{\alpha} in the posterior distribution of a
#' @param d_0_a parameter \eqn{d_0^a} in the posterior distribution of a
#' @param d_1_a parameter \eqn{d_1^a} in the posterior distribution of a
#' @param mu_star_n_r vector with number of observations allocated to each cluster
#' @param n.burn number of iterations in the simulation considered in the burn-in period.
#' @param n.thin number of iterations discarded between two simulated values (for thinning of the MCMC chain).
#' @param max.time maximum allowed time for the simulation process. The function returns \code{Error} if exceeded.
#' @param verbose if \code{T}, the function reports extra information on progress.
#' @param USING_CPP indicates usage of C++ in some modules.
#' 
#' @return A list with two elements:
#' \describe{
#'   \item{$a.chain}{A numeric vector with the simulated values from the posterior distribution of \emph{a}}
#'   \item{$accept.indic}{A numeric vector with the simulated values from the posterior distribution of \emph{a}}
#' }
#'
#' @references
#' Carmona C., Nieto-Barajas L., Canale A. (2017). \emph{Model based approach for household clustering with mixed scale variables.}
#'
#' @importFrom stats rbinom
#' @importFrom compiler cmpfun
#' 
#' @keywords internal

sampling_a <- function( n=1, a.ini,
                        b, alpha, d_0_a, d_1_a,
                        mu_star_n_r,
                        n.burn=0,n.thin=0,
                        max.time=Inf,
                        verbose=F,
                        USING_CPP=TRUE ) {
  
  # Metropolis-Hastings for 'a' #
  # MH Sampling from 'a' #
  
  if( a.ini<0 | a.ini>1 ){
    cat('\nError: the value for "a.ini" has to be in [0,1) \n')
    stop('the value for "a.ini" has to be in [0,1)') }
  if( (b+a.ini)<0 ){
    cat('\nError: the value for "b" has to be greater than -a.ini \n')
    stop('the value for "b" has to be greater than -a.ini')
  }
  
  # browser()
  
  #n.obs <- sum(mu_star_n_r)
  #r <- length(mu_star_n_r)
  
  # initializing the chain #
  a.chain <- as.numeric(NULL)
  accept.indic <- as.logical(NULL)
  a <- a.ini
  
  it_t_0 <- Sys.time()
  it_t_i <- as.numeric(Sys.time()-it_t_0)
  
  # total number of iterations
  n.iter <- n.burn + 1 + (n-1)*(n.thin+1)
  
  while( ( length(a.chain) < n.iter ) & (it_t_i<max.time) ) {
    
    if(verbose) {cat(".")}
    it_t_i <- as.numeric(Sys.time()-it_t_0)
    
    # generate proposal "a_new"
    if(b>0){
      a.prop <- rbinom( n=1, size=1, prob=1-alpha )
      a.prop <- ifelse( a.prop==0, 0, runif(n=1,min=0,max=1) )
    } else {
      a.prop <- runif(n=1,min=max(0,-b),max=1)
    }
    
    
    # posterior probability for the current value of a
    if(USING_CPP) {
      log_f_post_a <- log_f_post_a_cpp( a=a,
                                        b=b,
                                        alpha=alpha,
                                        d_0_a=d_0_a,d_1_a=d_1_a,
                                        mu_star_n_r=mu_star_n_r )
      
      log_f_post_a.prop <- log_f_post_a_cpp( a=a.prop,
                                             b=b,
                                             alpha=alpha,
                                             d_0_a=d_0_a,d_1_a=d_1_a,
                                             mu_star_n_r=mu_star_n_r )
    } else {
      log_f_post_a <- log_f_post_a( a=a,
                                    b=b,
                                    alpha=alpha,
                                    d_0_a=d_0_a,d_1_a=d_1_a,
                                    mu_star_n_r=mu_star_n_r )
      
      log_f_post_a.prop <- log_f_post_a( a=a.prop,
                                         b=b,
                                         alpha=alpha,
                                         d_0_a=d_0_a,d_1_a=d_1_a,
                                         mu_star_n_r=mu_star_n_r )
    }
    # posterior probability for the proposal value of a
    
    
    log_r <- log_f_post_a.prop - log_f_post_a
    # exp(log_r)
    
    #if( runif(1,0,1) > lik_ratio ) {
    if(is.na(log_r)) {browser()}
    if( log_r>=0 || runif(1,0,1) < exp(log_r) ) {
      a <- a.prop
      a.chain <- c( a.chain, a.prop )
      accept.indic <- c(accept.indic,T)
    } else {
      a.chain <- c( a.chain, a )
      accept.indic <- c(accept.indic,F)
    }
  }
  
  if(it_t_i>max.time) {
    cat('Error: TIMEOUT, There is a problem simulating from "a" \n')
    stop('TIMEOUT, There is a problem simulating from "a"')
  }
  
  index.iter.out <- seq(from=n.burn+1,to=n.iter,by=(n.thin+1))
  
  a.sim <- list( a.chain = a.chain[index.iter.out],
                 accept.indic = accept.indic )
  
  return( a.sim )
}

sampling_a <- compiler::cmpfun(sampling_a)