R/partitionSampler.R

In eifer4/stochasticSampling: Stochastic partitioning sampler.

partitionSampler <- function(componentData, priors, invariantData, n.iter, total.iter, iterCond, temp = 1, tempIdx) {
  
  iter.range <- seq(total.iter + 1, n.iter+total.iter, length.out = n.iter)
  progress.iter <- iterCond$progress.iter
  save.iter <- iterCond$save.iter
  n.burnin <- iterCond$n.burnin
  half.switch <- iterCond$half.switch
  n.save <- ceiling(sum(iter.range>n.burnin)/save.iter)
  if(tempIdx ==1){
    save <- TRUE
  } else {
    save <- FALSE
  }

  Q <- invariantData$Q
  P <- invariantData$P
  Y_corr  <- array(0, dim=c(Q, Q, 2), dimnames=list(yrows=NULL, ycols=NULL, half=1:2))
  X_corr  <- array(0, dim=c(Q, P, 2), dimnames=list(yrows=NULL, xcols=NULL, half=1:2))
  log_prob <- numComp <-  numeric(n.save)
  
  accept <- numeric(2)
  hidx <- sidx <- 1
  
  out <- componentData
  
  for(iter in iter.range) {

    out <- updateParam(comp=out$comp, priors=priors,  invariantData=invariantData, compCount=out$compCount, temp=temp)
    
    accept <- accept + out$accept
    
    if(iter > n.burnin & save){
      for(k in 1:length(out$comp)) {
        Y_corr[out$comp[[k]]$Y, out$comp[[k]]$Y, hidx] <- Y_corr[out$comp[[k]]$Y,out$comp[[k]]$Y, hidx] + 1
        X_corr[out$comp[[k]]$Y, out$comp[[k]]$X, hidx] <- X_corr[out$comp[[k]]$Y,out$comp[[k]]$X, hidx] + 1
      }
      if(iter == half.switch) hidx <- hidx + 1
      
      if((iter %% save.iter)==0) {
        
        log_prob[sidx] <- logProb(out$comp, invariantData, priors, temp=1.0)
        numComp[sidx] <- out$compCount$K
        
        sidx <- sidx + 1
      }
    }
    
  }
  
  return(
    list(Y_corr = Y_corr,
         X_corr = X_corr,
         log_prob = log_prob,
         numComp = numComp,
         componentData = list(comp = out$comp, compCount = out$compCount),
         n.save=n.save,
         accept = accept)
  )
}