R/quantile.residuals.DDPstar.R

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

quantileResiduals.DDPstar <- 
function(object, parallel = c("no", "multicore", "snow"), ncpus = 1, cl = NULL,...) {

	doMCMC <- function(k, object) {
		ncomp <- ncol(object$fit$probs)
		cdf <- rep(0, length(object$data_model$y))
		if(ncomp == 1) {
			mu <- object$data_model$X%*%matrix(object$fit$beta[k,,], ncol = 1)
		} else {
			mu <- object$data_model$X%*%t(object$fit$beta[k,,])
		}
		w <- object$fit$probs[k,]
		sd <- object$fit$sd[k,]
		for(i in 1:length(object$data_model$y)){
	  		cdf[i] <- sum(w*pnorm(object$data_model$y[i], mu[i,], sd))
		}
		res <- quantile(qnorm(cdf), ppoints(length(object$data_model$y)), type = 2)
		res
	}
	if(class(object)[1] != "DDPstar") {
		stop(paste0("This function can not be used for this object class: ", class(object)[object]))
	}

	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, 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)
                        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)
                        }                        
                    }
                }
            } else {
                lapply(seq_len(nsim), doMCMC, object = object)
            }

        quantile_residuals <- simplify2array(resMCMC)

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

	trajhat_m  <- apply(quantile_residuals, 1, mean)
	trajhat_ql <- apply(quantile_residuals, 1, quantile, prob = 0.025)
	trajhat_qh <- apply(quantile_residuals, 1, quantile, prob = 0.975)

	grid <- qnorm(ppoints(length(object$data_model$y)))
	df_quant_res <- data.frame(x = grid, y = trajhat_m, ymin = trajhat_ql, ymax = trajhat_qh)

	res <- list()
	res$x <- grid
	res$y <- trajhat_m
	res$ymin <- trajhat_ql
	res$ymax <- trajhat_qh
	invisible(res)
}