Nothing
R/RcppExports.R
In rjpdmp: Reversible Jump PDMP Samplers
Defines functions
zigzag_rr
zigzag_logit_ss
zigzag_logit
gibbs_logit
bps_n_rr
bps_s_rr
bps_s_logit
bps_n_logit
Documented in
bps_n_logit
bps_n_rr
bps_s_logit
bps_s_rr
gibbs_logit
zigzag_logit
zigzag_logit_ss
zigzag_rr
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
#' bps_n_logit
#'
#' Applies the Reversible Jump BPS Sampler with Velocities distributed from the Normal distribution to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler (Default has 1s on all components). This should be chosen with unit velocities on each component (regardless of sign).
#' @param ref Double for the refreshment rate of the BPS.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'\dontrun{
#' bps_fit <- bps_n_logit(maxTime = 1, dataX = data$dataX, datay = data$dataY,
#' prior_sigma2 = 10, theta0 = rep(0, p),
#' x0 = rep(0, p), ref = 0.1, rj_val = 0.6,
#' ppi = ppi, nmax = 1e6, burn = -1)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay =data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma =rep(0,p),
#' ppi = ppi)
#'
#' plot_pdmp(bps_fit, coords = 1:2, inds = 1:1e3,burn = .1, nsamples = 1e4,
#' mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
bps_n_logit <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, ref = 0.1, rj_val = 0.6, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_bps_n_logit`, maxTime, dataX, datay, prior_sigma2, x0, theta0, ref, rj_val, ppi, nmax, burn)
}
#' bps_s_logit
#'
#' Applies the Reversible Jump BPS Sampler with Velocities drawn Uniformly on the p-Sphere to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler. This should be chosen with unit velocities on each component (regardless of sign).
#' @param ref Double for the refreshment rate of the BPS.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'
#'\dontrun{
#' bps_fit <- bps_s_logit(maxTime = 1, dataX = data$dataX, datay = data$dataY,
#' prior_sigma2 = 10, theta0 = rep(0, p),
#' x0 = rep(0, p), ref = 0.1, rj_val = 0.6,
#' ppi = ppi)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay =data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma =rep(0,p),
#' ppi = ppi)
#'
#' plot_pdmp(bps_fit, coords = 1:2, inds = 1:1e4,burn = .1, nsamples = 1e4,
#' mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
bps_s_logit <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, ref = 0.01, rj_val = 0.6, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_bps_s_logit`, maxTime, dataX, datay, prior_sigma2, x0, theta0, ref, rj_val, ppi, nmax, burn)
}
#' bps_s_rr
#'
#' Applies the Reversible Jump BPS Sampler with Velocities drawn Uniformly on the p-Sphere to a Robust Regression problem with dirac spike and slab prior.
#' Included variables are given an independent Gaussian prior with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a continuous response variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler.
#' @param ref Refreshment rate for BPS.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burn-in. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switches velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#'
#' generate.rr.data <- function(beta, n, Sig, noise, interc = TRUE) {
#' p <- length(beta)-(interc == TRUE)
#' dataX <- MASS::mvrnorm(n=n,mu=rep(0,p),Sigma=Sig)
#' if(interc) {dataX <- cbind(1, dataX)}
#' dataY <- rep(0, n)
#' dataY <- dataX %*% as.vector(beta)+rnorm(n, sd = sqrt(noise))
#' return(list(dataX = dataX, dataY = dataY))
#' }
#' p <- 3;
#' n<- 120
#' beta <- c(0.5,0.5, rep(0,p-1))
#' set.seed(1)
#' data <- generate.rr.data(beta,n,diag(1,p+1), noise = 2, interc = FALSE)
#' dataX <- data$dataX; dataY <- data$dataY
#'\dontrun{
#' set.seed(1)
#' ppi_val <- 1/4
#' res <- bps_s_rr(maxTime = 1, dataX = dataX, datay = dataY,
#' prior_sigma2 = 1e2, x0 = rep(0,p+1), theta0 = rep(0,p+1),
#' rj_val = 0.6, ppi = ppi_val, nmax = 1e5)
#'
#' plot_pdmp(res, coords = 1:3, inds = 1:1e3)
#'}
#' @export
bps_s_rr <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, ref = 0.1, rj_val = 0.5, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_bps_s_rr`, maxTime, dataX, datay, prior_sigma2, x0, theta0, ref, rj_val, ppi, nmax, burn)
}
#' bps_n_rr
#'
#' Applies the Reversible Jump BPS Sampler with Velocities drawn from the Normal distribution to a Robust Regression problem with dirac spike and slab prior.
#' Included variables are given an independent Gaussian prior with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a continuous response variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler.
#' @param ref Refreshment rate for BPS.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.rr.data <- function(beta, n, Sig, noise, interc = TRUE) {
#' p <- length(beta)-(interc == TRUE)
#' dataX <- MASS::mvrnorm(n=n,mu=rep(0,p),Sigma=Sig)
#' if(interc) {dataX <- cbind(1, dataX)}
#' dataY <- rep(0, n)
#' dataY <- dataX %*% as.vector(beta)+rnorm(n, sd = sqrt(noise))
#' return(list(dataX = dataX, dataY = dataY))
#' }
#' p <- 3;
#' n<- 120
#' beta <- c(0.5,0.5, rep(0,p-1))
#' set.seed(1)
#' data <- generate.rr.data(beta,n,diag(1,p+1), noise = 2, interc = FALSE)
#' dataX <- data$dataX; dataY <- data$dataY
#'\dontrun{
#' set.seed(1)
#' ppi_val <- 1/4
#' res <- bps_n_rr(maxTime = 1, dataX = dataX, datay = dataY,
#' prior_sigma2 = 1e2, x0 = rep(0,p+1), theta0 = rep(0,p+1),
#' rj_val = 0.6, ppi = ppi_val, nmax = 1e5, ref = 0.1, burn = -1)
#'
#' plot_pdmp(res, coords = 1:3, inds = 1:1e3)
#'}
#' @export
bps_n_rr <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, ref = 0.1, rj_val = 0.5, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_bps_n_rr`, maxTime, dataX, datay, prior_sigma2, x0, theta0, ref, rj_val, ppi, nmax, burn)
}
#' gibbs_logit
#'
#' Applies the Collapsed Gibbs Sampler to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#' Code makes use of the package set-up for Polya-Gamma simulation available at \code{https://github.com/jgscott/helloPG}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables. Default 10.
#' @param beta Initial position of the regression parameter
#' @param gamma Initial model for the sampler. Enteries should either be 1s or 0s.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nsamples Maximum number of samples. Default value is 10^5, lower values should be chosen for memory constraints if less samples are desired.
#' @return Returns a list with the following objects:
#' @return \code{beta}: Matrix of regression parameter samples, columns are samples.
#' @return \code{gamma}: Matrix of model parameter samples columns are samples.
#' @return \code{times}: computation times at sampled events - Useful for plotting computational efficiency.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'
#' zigzag_fit <- zigzag_logit(maxTime = 1, dataX = data$dataX,
#' datay = data$dataY, prior_sigma2 = 10,
#' theta0 = rep(0, p), x0 = rep(0, p), rj_val = 0.6,
#' ppi = ppi)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay =data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma =rep(0,p),
#' ppi = ppi)
#'\dontrun{
#' plot_pdmp(zigzag_fit, coords = 1:2, inds = 1:1e3,burn = .1,
#' nsamples = 1e4, mcmc_samples =t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
gibbs_logit <- function(dataX, datay, beta, gamma, ppi = 0.5, nsamples = 1e5L, maxTime = 1e8, prior_sigma2 = 10.0) {
.Call(`_rjpdmp_gibbs_logit`, dataX, datay, beta, gamma, ppi, nsamples, maxTime, prior_sigma2)
}
#' zigzag_logit
#'
#' Applies the Reversible Jump ZigZag Sampler to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler (Default has 1s on all components). This should be chosen with unit velocities on each component (regardless of sign).
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'
#' zigzag_fit <- zigzag_logit(maxTime = 1, dataX = data$dataX, datay = data$dataY,
#' prior_sigma2 = 10,theta0 = rep(0, p), x0 = rep(0, p), rj_val = 0.6,
#' ppi = ppi)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay = data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma = rep(0,p),
#' ppi = ppi)
#'\dontrun{
#' plot_pdmp(zigzag_fit, coords = 1:2, inds = 1:1e3,burn = .1, nsamples = 1e4,
#' mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
zigzag_logit <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, rj_val = 0.6, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_zigzag_logit`, maxTime, dataX, datay, prior_sigma2, x0, theta0, rj_val, ppi, nmax, burn)
}
#' zigzag_logit_ss
#'
#' Applies the Reversible Jump ZigZag Sampler with subsampling to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler (Default has 1s on all components). This should be chosen with unit velocities on each component (regardless of sign).
#' @param cvref Control variate vector of dimension p for subsampling. If no control variate set to a vector of zeros.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'
#'\dontrun{
#' zigzag_fit <- zigzag_logit(maxTime = 1, dataX = data$dataX,
#' datay = data$dataY, prior_sigma2 = 10,
#' theta0 = rep(0, p), x0 = rep(0, p), rj_val = 0.6,
#' ppi = ppi)
#'
#' zigzag_fit_s <- zigzag_logit_ss(maxTime = 1, dataX = data$dataX,
#' datay = data$dataY,prior_sigma2 = 10,
#' theta0 = rep(0, p), x0 = rep(0, p),
#' rj_val = 0.6, cvref = c(1,rep(0,p-1)),
#' ppi = ppi)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay =data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma =rep(0,p),
#' ppi = ppi)
#'
#' plot_pdmp_multiple(list(zigzag_fit,zigzag_fit_s), coords = 1:2, burn = .1,
#' inds = 1:1e2, nsamples = 1e4,
#' mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
zigzag_logit_ss <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, cvref, rj_val = 0.6, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_zigzag_logit_ss`, maxTime, dataX, datay, prior_sigma2, x0, theta0, cvref, rj_val, ppi, nmax, burn)
}
#' zigzag_rr
#'
#' Applies the Reversible Jump ZigZag Sampler to a Robust Regression problem with dirac spike and slab prior.
#' Included variables are given an independent Gaussian prior with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a continuous response variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler (Default has 1s on all components). This should be chosen with unit velocities on each component (regardless of sign).
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.rr.data <- function(beta, n, Sig, noise, interc = TRUE) {
#' p <- length(beta)-(interc == TRUE)
#' dataX <- MASS::mvrnorm(n=n,mu=rep(0,p),Sigma=Sig)
#' if(interc) {dataX <- cbind(1, dataX)}
#' dataY <- rep(0, n)
#' dataY <- dataX %*% as.vector(beta)+rnorm(n, sd = sqrt(noise))
#' return(list(dataX = dataX, dataY = dataY))
#' }
#' p <- 3;
#' n<- 120
#' beta <- c(0.5,0.5, rep(0,p-1))
#' set.seed(1)
#' data <- generate.rr.data(beta,n,diag(1,p+1), noise = 2, interc = FALSE)
#' dataX <- data$dataX; dataY <- data$dataY
#'
#' set.seed(1)
#' ppi_val <- 1/4
#' res <- zigzag_rr(maxTime = 1, dataX = dataX, datay = dataY,
#' prior_sigma2 = 1e2, x0 = rep(0,p+1), theta0 = rep(0,p+1),
#' rj_val = 0.6, ppi = ppi_val, nmax = 1e5)
#'\dontrun{
#' plot_pdmp(res, coords = 1:3, inds = 1:1e3)
#'}
#'
#' @export
zigzag_rr <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, rj_val = 0.5, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_zigzag_rr`, maxTime, dataX, datay, prior_sigma2, x0, theta0, rj_val, ppi, nmax, burn)
}
Try the rjpdmp package in your browser
Any scripts or data that you put into this service are public.
rjpdmp documentation built on March 18, 2022, 7:52 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions zigzag_rr zigzag_logit_ss zigzag_logit gibbs_logit bps_n_rr bps_s_rr bps_s_logit bps_n_logit
Documented in bps_n_logit bps_n_rr bps_s_logit bps_s_rr gibbs_logit zigzag_logit zigzag_logit_ss zigzag_rr
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
#' bps_n_logit
#'
#' Applies the Reversible Jump BPS Sampler with Velocities distributed from the Normal distribution to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler (Default has 1s on all components). This should be chosen with unit velocities on each component (regardless of sign).
#' @param ref Double for the refreshment rate of the BPS.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'\dontrun{
#' bps_fit <- bps_n_logit(maxTime = 1, dataX = data$dataX, datay = data$dataY,
#' prior_sigma2 = 10, theta0 = rep(0, p),
#' x0 = rep(0, p), ref = 0.1, rj_val = 0.6,
#' ppi = ppi, nmax = 1e6, burn = -1)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay =data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma =rep(0,p),
#' ppi = ppi)
#'
#' plot_pdmp(bps_fit, coords = 1:2, inds = 1:1e3,burn = .1, nsamples = 1e4,
#' mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
bps_n_logit <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, ref = 0.1, rj_val = 0.6, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_bps_n_logit`, maxTime, dataX, datay, prior_sigma2, x0, theta0, ref, rj_val, ppi, nmax, burn)
}
#' bps_s_logit
#'
#' Applies the Reversible Jump BPS Sampler with Velocities drawn Uniformly on the p-Sphere to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler. This should be chosen with unit velocities on each component (regardless of sign).
#' @param ref Double for the refreshment rate of the BPS.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'
#'\dontrun{
#' bps_fit <- bps_s_logit(maxTime = 1, dataX = data$dataX, datay = data$dataY,
#' prior_sigma2 = 10, theta0 = rep(0, p),
#' x0 = rep(0, p), ref = 0.1, rj_val = 0.6,
#' ppi = ppi)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay =data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma =rep(0,p),
#' ppi = ppi)
#'
#' plot_pdmp(bps_fit, coords = 1:2, inds = 1:1e4,burn = .1, nsamples = 1e4,
#' mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
bps_s_logit <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, ref = 0.01, rj_val = 0.6, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_bps_s_logit`, maxTime, dataX, datay, prior_sigma2, x0, theta0, ref, rj_val, ppi, nmax, burn)
}
#' bps_s_rr
#'
#' Applies the Reversible Jump BPS Sampler with Velocities drawn Uniformly on the p-Sphere to a Robust Regression problem with dirac spike and slab prior.
#' Included variables are given an independent Gaussian prior with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a continuous response variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler.
#' @param ref Refreshment rate for BPS.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burn-in. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switches velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#'
#' generate.rr.data <- function(beta, n, Sig, noise, interc = TRUE) {
#' p <- length(beta)-(interc == TRUE)
#' dataX <- MASS::mvrnorm(n=n,mu=rep(0,p),Sigma=Sig)
#' if(interc) {dataX <- cbind(1, dataX)}
#' dataY <- rep(0, n)
#' dataY <- dataX %*% as.vector(beta)+rnorm(n, sd = sqrt(noise))
#' return(list(dataX = dataX, dataY = dataY))
#' }
#' p <- 3;
#' n<- 120
#' beta <- c(0.5,0.5, rep(0,p-1))
#' set.seed(1)
#' data <- generate.rr.data(beta,n,diag(1,p+1), noise = 2, interc = FALSE)
#' dataX <- data$dataX; dataY <- data$dataY
#'\dontrun{
#' set.seed(1)
#' ppi_val <- 1/4
#' res <- bps_s_rr(maxTime = 1, dataX = dataX, datay = dataY,
#' prior_sigma2 = 1e2, x0 = rep(0,p+1), theta0 = rep(0,p+1),
#' rj_val = 0.6, ppi = ppi_val, nmax = 1e5)
#'
#' plot_pdmp(res, coords = 1:3, inds = 1:1e3)
#'}
#' @export
bps_s_rr <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, ref = 0.1, rj_val = 0.5, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_bps_s_rr`, maxTime, dataX, datay, prior_sigma2, x0, theta0, ref, rj_val, ppi, nmax, burn)
}
#' bps_n_rr
#'
#' Applies the Reversible Jump BPS Sampler with Velocities drawn from the Normal distribution to a Robust Regression problem with dirac spike and slab prior.
#' Included variables are given an independent Gaussian prior with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a continuous response variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler.
#' @param ref Refreshment rate for BPS.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.rr.data <- function(beta, n, Sig, noise, interc = TRUE) {
#' p <- length(beta)-(interc == TRUE)
#' dataX <- MASS::mvrnorm(n=n,mu=rep(0,p),Sigma=Sig)
#' if(interc) {dataX <- cbind(1, dataX)}
#' dataY <- rep(0, n)
#' dataY <- dataX %*% as.vector(beta)+rnorm(n, sd = sqrt(noise))
#' return(list(dataX = dataX, dataY = dataY))
#' }
#' p <- 3;
#' n<- 120
#' beta <- c(0.5,0.5, rep(0,p-1))
#' set.seed(1)
#' data <- generate.rr.data(beta,n,diag(1,p+1), noise = 2, interc = FALSE)
#' dataX <- data$dataX; dataY <- data$dataY
#'\dontrun{
#' set.seed(1)
#' ppi_val <- 1/4
#' res <- bps_n_rr(maxTime = 1, dataX = dataX, datay = dataY,
#' prior_sigma2 = 1e2, x0 = rep(0,p+1), theta0 = rep(0,p+1),
#' rj_val = 0.6, ppi = ppi_val, nmax = 1e5, ref = 0.1, burn = -1)
#'
#' plot_pdmp(res, coords = 1:3, inds = 1:1e3)
#'}
#' @export
bps_n_rr <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, ref = 0.1, rj_val = 0.5, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_bps_n_rr`, maxTime, dataX, datay, prior_sigma2, x0, theta0, ref, rj_val, ppi, nmax, burn)
}
#' gibbs_logit
#'
#' Applies the Collapsed Gibbs Sampler to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#' Code makes use of the package set-up for Polya-Gamma simulation available at \code{https://github.com/jgscott/helloPG}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables. Default 10.
#' @param beta Initial position of the regression parameter
#' @param gamma Initial model for the sampler. Enteries should either be 1s or 0s.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nsamples Maximum number of samples. Default value is 10^5, lower values should be chosen for memory constraints if less samples are desired.
#' @return Returns a list with the following objects:
#' @return \code{beta}: Matrix of regression parameter samples, columns are samples.
#' @return \code{gamma}: Matrix of model parameter samples columns are samples.
#' @return \code{times}: computation times at sampled events - Useful for plotting computational efficiency.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'
#' zigzag_fit <- zigzag_logit(maxTime = 1, dataX = data$dataX,
#' datay = data$dataY, prior_sigma2 = 10,
#' theta0 = rep(0, p), x0 = rep(0, p), rj_val = 0.6,
#' ppi = ppi)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay =data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma =rep(0,p),
#' ppi = ppi)
#'\dontrun{
#' plot_pdmp(zigzag_fit, coords = 1:2, inds = 1:1e3,burn = .1,
#' nsamples = 1e4, mcmc_samples =t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
gibbs_logit <- function(dataX, datay, beta, gamma, ppi = 0.5, nsamples = 1e5L, maxTime = 1e8, prior_sigma2 = 10.0) {
.Call(`_rjpdmp_gibbs_logit`, dataX, datay, beta, gamma, ppi, nsamples, maxTime, prior_sigma2)
}
#' zigzag_logit
#'
#' Applies the Reversible Jump ZigZag Sampler to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler (Default has 1s on all components). This should be chosen with unit velocities on each component (regardless of sign).
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'
#' zigzag_fit <- zigzag_logit(maxTime = 1, dataX = data$dataX, datay = data$dataY,
#' prior_sigma2 = 10,theta0 = rep(0, p), x0 = rep(0, p), rj_val = 0.6,
#' ppi = ppi)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay = data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma = rep(0,p),
#' ppi = ppi)
#'\dontrun{
#' plot_pdmp(zigzag_fit, coords = 1:2, inds = 1:1e3,burn = .1, nsamples = 1e4,
#' mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
zigzag_logit <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, rj_val = 0.6, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_zigzag_logit`, maxTime, dataX, datay, prior_sigma2, x0, theta0, rj_val, ppi, nmax, burn)
}
#' zigzag_logit_ss
#'
#' Applies the Reversible Jump ZigZag Sampler with subsampling to a Logistic regression with dirac spike and slab distribution, as detailed in Reversible Jump PDMP Samplers for Variable Selection, 2020.
#' For included variables an independent Gaussian prior is assumed with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a {0, 1}-valued variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler (Default has 1s on all components). This should be chosen with unit velocities on each component (regardless of sign).
#' @param cvref Control variate vector of dimension p for subsampling. If no control variate set to a vector of zeros.
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.logistic.data <- function(beta, n.obs, Sig) {
#' p <- length(beta)
#' dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
#' vals <- dataX %*% as.vector(beta)
#' generateY <- function(p) { rbinom(1, 1, p)}
#' dataY <- sapply(1/(1 + exp(-vals)), generateY)
#' return(list(dataX = dataX, dataY = dataY))
#' }
#'
#' n <- 15
#' p <- 25
#' beta <- c(1, rep(0, p-1))
#' Siginv <- diag(1,p,p)
#' Siginv[1,2] <- Siginv[2,1] <- 0.9
#' set.seed(1)
#' data <- generate.logistic.data(beta, n, solve(Siginv))
#' ppi <- 2/p
#'
#'\dontrun{
#' zigzag_fit <- zigzag_logit(maxTime = 1, dataX = data$dataX,
#' datay = data$dataY, prior_sigma2 = 10,
#' theta0 = rep(0, p), x0 = rep(0, p), rj_val = 0.6,
#' ppi = ppi)
#'
#' zigzag_fit_s <- zigzag_logit_ss(maxTime = 1, dataX = data$dataX,
#' datay = data$dataY,prior_sigma2 = 10,
#' theta0 = rep(0, p), x0 = rep(0, p),
#' rj_val = 0.6, cvref = c(1,rep(0,p-1)),
#' ppi = ppi)
#'
#' gibbs_fit <- gibbs_logit(maxTime = 1, dataX = data$dataX, datay =data$dataY,
#' prior_sigma2 = 10,beta = rep(0,p), gamma =rep(0,p),
#' ppi = ppi)
#'
#' plot_pdmp_multiple(list(zigzag_fit,zigzag_fit_s), coords = 1:2, burn = .1,
#' inds = 1:1e2, nsamples = 1e4,
#' mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))
#'}
#' @export
zigzag_logit_ss <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, cvref, rj_val = 0.6, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_zigzag_logit_ss`, maxTime, dataX, datay, prior_sigma2, x0, theta0, cvref, rj_val, ppi, nmax, burn)
}
#' zigzag_rr
#'
#' Applies the Reversible Jump ZigZag Sampler to a Robust Regression problem with dirac spike and slab prior.
#' Included variables are given an independent Gaussian prior with variance \code{prior_sigma2} and mean zero, variables are given prior probabilities of inclusion \code{ppi}.
#'
#' @param maxTime Maximum runtime (in Seconds) of the algorithm; will terminate the code after a given computation time or nmax iterations of the algorithm is reached.
#' @param dataX Matrix of all covariates where the i-th row corresponds to all p covariates x_{i,1}, ..., x_{i,p} of the i-th observation.
#' @param datay Vector of n observations of a continuous response variable y.
#' @param prior_sigma2 Double for the prior variance for included variables.
#' @param x0 Initial position of the regression parameter
#' @param theta0 Initial velocity for the sampler (Default has 1s on all components). This should be chosen with unit velocities on each component (regardless of sign).
#' @param rj_val Reversible jump parameter for the PDMP method. This value is fixed over all models and is interpreted as the probability to jump to a reduced model when a parameter hits zero.
#' @param ppi Double for the prior probability of inclusion (ppi) for each parameter.
#' @param nmax Maximum number of iterations (simulated events) of the algorithm; will stop the algorithm when this number of iterations of the method have occured. Default value is 1e6, lower values should be chosen for memory constraints if less iterations are desired.
#' @param burn Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.
#' @return Returns a list with the following objects:
#' @return \code{times}: Vector of event times where ZigZag process switchs velocity or jumps models.
#' @return \code{positions}: Matrix of positions at which event times occur, these are not samples from the PDMP.
#' @return \code{theta}: Matrix of new velocities at event times.
#' @examples
#'
#' generate.rr.data <- function(beta, n, Sig, noise, interc = TRUE) {
#' p <- length(beta)-(interc == TRUE)
#' dataX <- MASS::mvrnorm(n=n,mu=rep(0,p),Sigma=Sig)
#' if(interc) {dataX <- cbind(1, dataX)}
#' dataY <- rep(0, n)
#' dataY <- dataX %*% as.vector(beta)+rnorm(n, sd = sqrt(noise))
#' return(list(dataX = dataX, dataY = dataY))
#' }
#' p <- 3;
#' n<- 120
#' beta <- c(0.5,0.5, rep(0,p-1))
#' set.seed(1)
#' data <- generate.rr.data(beta,n,diag(1,p+1), noise = 2, interc = FALSE)
#' dataX <- data$dataX; dataY <- data$dataY
#'
#' set.seed(1)
#' ppi_val <- 1/4
#' res <- zigzag_rr(maxTime = 1, dataX = dataX, datay = dataY,
#' prior_sigma2 = 1e2, x0 = rep(0,p+1), theta0 = rep(0,p+1),
#' rj_val = 0.6, ppi = ppi_val, nmax = 1e5)
#'\dontrun{
#' plot_pdmp(res, coords = 1:3, inds = 1:1e3)
#'}
#'
#' @export
zigzag_rr <- function(maxTime, dataX, datay, prior_sigma2, x0, theta0, rj_val = 0.5, ppi = 0.5, nmax = 1e6L, burn = -1L) {
.Call(`_rjpdmp_zigzag_rr`, maxTime, dataX, datay, prior_sigma2, x0, theta0, rj_val, ppi, nmax, burn)
}
Try the rjpdmp package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.