R/predict.DDPstar.probfun.R

In DDPstar: Density Regression via Dirichlet Process Mixtures of Normal Structured Additive Regression Models

predict_DDPstar.probfun <- function(object, newdata, q.value, parallel = c("no", "multicore", "snow"), ncpus = 1, cl = NULL) {
  
  doMCMC <- function(k, object, Xp, q.value) {
    probfun <- matrix(0, nrow = length(q.value), ncol = nrow(Xp))
    ncomp <- ncol(object$fit$probs)
    if(ncomp == 1) {
      Beta <- matrix(object$fit$beta[k,,], nrow = 1)
    } else {
      Beta <- object$fit$beta[k,,]
    }

    mu <- Xp%*%t(Beta)
    sigma <- object$fit$sd[k,]
    w <- object$fit$probs[k,]
    
    for(i in 1:length(q.value)) {
      for(j in 1:nrow(Xp)) {
        probfun[i,j] <- sum(w*pnorm(q.value[i], mean = c(mu[j,]), sd = sigma))
      }
    }
    res <- list()
    res$probfun <- probfun
    res$foo <- 1
    res
  }

  Xp <- predict.design.matrix.DDPstar(object = object, newdata = newdata)$X
  
  parallel <- match.arg(parallel)

  nsim <- nrow(object$fit$sd)
  if(nsim > 0) {
    do_mc <- do_snow <- FALSE
      if (parallel != "no" && ncpus > 1L) {
          if (parallel == "multicore") {
              do_mc <- .Platform$OS.type != "windows"
          } else if (parallel == "snow") {
              do_snow <- TRUE
          }
          if (!do_mc && !do_snow) {
              ncpus <- 1L
          }       
          loadNamespace("parallel") # get this out of the way before recording seed
      }

      resMCMC <- if (ncpus > 1L && (do_mc || do_snow)) {
                if (do_mc) {
                    parallel::mclapply(seq_len(nsim), doMCMC, object = object, Xp = Xp, q.value = q.value, mc.cores = ncpus)
                } else if (do_snow) {                
                    if (is.null(cl)) {
                        cl <- parallel::makePSOCKcluster(rep("localhost", ncpus))
                        if(RNGkind()[1L] == "L'Ecuyer-CMRG") {
                            parallel::clusterSetRNGStream(cl)
                        }
                        res <- parallel::parLapply(cl, seq_len(nsim), doMCMC, object = object, Xp = Xp, q.value = q.value)
                        parallel::stopCluster(cl)
                        res
                    } else {
                        if(!inherits(cl, "cluster")) {
                            stop("Class of object 'cl' is not correct")
                        } else {
                            parallel::parLapply(cl, seq_len(nsim), doMCMC, object = object, Xp = Xp, q.value = q.value)
                        }                        
                    }
                }
            } else {
                lapply(seq_len(nsim), doMCMC, object = object, Xp = Xp, q.value = q.value)
            }

        resMCMC <- simplify2array(resMCMC)
        probfun<- simplify2array(resMCMC["probfun",])

  } else {
    stop("nsave should be larger than zero.")
  }

  res <- list()
  res$Xp <- Xp
  res$q.value <- q.value
  res$probfun <- probfun
  class(res) <- "DDPstar.probfun"
  res
}