gen_sample: Generate samples for PDMP trajectory

In rjpdmp: Reversible Jump PDMP Samplers

Generate samples for PDMP trajectory

Description

Get samples from PDMP trajectories taking a fixed time discretisation.

Usage

gen_sample(positions, times, nsample, theta = NULL, burn = 1)

Arguments

positions	Matrix of positions from the PDMP trajectory, each column should correspond to a position
times	Vector of event times from the PDMP trajectory
nsample	Number of desired samples from the PDMP trajectory
theta	Optional Matrix of velocities from the PDMP trajectory, each column should correspond to a velocity
burn	Index to start the discretisation from. Default is 1.

Value

Returns a list with the following objects:

x: Matrix of extracted samples of the position (x) taken using a fixed time discretisation of the PDMP

theta: Matrix of extracted samples of the velocity (theta) taken using a fixed time discretisation of the PDMP

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: 
samples <- gen_sample(zigzag_fit$positions, zigzag_fit$times, 10^4)

plot(zigzag_fit$positions[1,],zigzag_fit$positions[2,], type = 'l', xlab = 'x1', ylab = 'x2')
points(samples$xx[1,], samples$xx[2,], col='red', pch=20)

## End(Not run)

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