rjpdmp (version 2.0.0)

bps_s_logit: bps_s_logit

Description

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 prior_sigma2 and mean zero, variables are given prior probabilities of inclusion ppi.

Usage

bps_s_logit(
  maxTime,
  dataX,
  datay,
  prior_sigma2,
  x0,
  theta0,
  ref = 0.01,
  rj_val = 0.6,
  ppi = 0.5,
  nmax = 1000000L,
  burn = -1L
)

Arguments

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.

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.

datay

Vector of n observations of a 0, 1-valued variable y.

prior_sigma2

Double for the prior variance for included variables.

x0

Initial position of the regression parameter

theta0

Initial velocity for the sampler. This should be chosen with unit velocities on each component (regardless of sign).

ref

Double for the refreshment rate of the BPS.

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.

ppi

Double for the prior probability of inclusion (ppi) for each parameter.

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.

burn

Optional number of iterations to use for burnin. These are not stored so can be useful in memory intensive problems.

Value

Returns a list with the following objects:

times: Vector of event times where ZigZag process switchs velocity or jumps models.

positions: Matrix of positions at which event times occur, these are not samples from the PDMP.

theta: Matrix of new velocities at event times.

Examples

# NOT RUN {
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

# }
# NOT RUN {
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))
# }