model_probabilities: Calculate posterior probabilities of inclusion based on PDMP...

In rjpdmp: Reversible Jump PDMP Samplers

Calculate posterior probabilities of inclusion based on PDMP trajectories

Description

Calculate either marginal probabilities of inclusions or posterior probabilities of specific models.

Usage

model_probabilities(times, thetas, models = NULL, marginals = NULL, burnin = 1)

Arguments

times	Vector of event times from the PDMP trajectory
thetas	Matrix of velocities from the PDMP trajectory, each column should correspond to a velocities
models	Optional Matrix of indicies where rows correspond to models. Will return probabilities of each model prob_mod.
marginals	Optional Vector of indices to calculate the marginal probabilities of inclusion. Will return probabilities of inclusion for variable index marginal_prob.
burnin	Number of events to use as burnin

Value

Returns a list with the following objects:

prob_mod: Vector of posterior model probabilities based on the PDMP trajectories

marginal_prob: Vector of marginal probabilities for inclusion

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)
## Not run: 
a <- models_visited(zigzag_fit$theta)

# Work out probability of top 10 most visited models and all marginal inclusion probabilities
# specific model probabilities become trivially small for large dimensions
b <- model_probabilities(zigzag_fit$times, zigzag_fit$thetas,
                         models = a[1:10,1:p], marginals=1:p)

## End(Not run)

rjpdmp documentation built on March 18, 2022, 7:52 p.m.