plot_pdmp: Plot PDMP dynamics and samples for posterior distributions

In rjpdmp: Reversible Jump PDMP Samplers

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

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))

## End(Not run)

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