rjpdmp (version 2.0.0)

plot_pdmp: Plot PDMP dynamics and samples for posterior distributions

Description

Plot marginal densities and joint pairs plots for trajectories and samples of PDMP samplers and optionally MCMC samples for comparison. Care should be taken when interpreting marginal KDE estimates on the diagonal as the bandwidth of the KDE has an impact on how the Dirac spike is visualised.

Usage

plot_pdmp(
  pdmp_res,
  coords = 1:2,
  inds = 1:10^3,
  nsamples = 10^3,
  burn = 0.1,
  mcmc_samples = NULL,
  pch = 20,
  cols = NULL
)

Arguments

pdmp_res

List of positions, times and velocities returned from a PDMP sampler

coords

Vector of coordinates to plot the marginal and joint distributions

inds

Vector of indices of the PDMP trajectories to plot.

nsamples

Number of samples to generate and use for marginal density estimates of the PDMP methods

burn

Percentage of events to use as burn-in. Should be between 0 and 1.

mcmc_samples

Optional Matrix of samples from an MCMC method. Each row should be a sample.

pch

The graphics parameter for off diagonal plots. Default is 20.

cols

Colours to be used for plotting the PDMPs and MCMC samples (in order).

Value

Generates a plot of the marginal density on the diagonal and pairs plots of the trajectories

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

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)
# }
# NOT RUN {
plot_pdmp(zigzag_fit, coords = 1:2, inds = 1:10^3,burn = .1,
          nsamples = 1e4, mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))
# }
# NOT RUN {
# }