| Title: | Zig-Zag Sampler |
|---|---|
| Description: | Implements the Zig-Zag algorithm (Bierkens, Fearnhead, Roberts, 2016) <arXiv:1607.03188> applied and Bouncy Particle Sampler <arXiv:1510.02451> for a Gaussian target and Student distribution. |
| Authors: | Joris Bierkens |
| Maintainer: | Joris Bierkens <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.2.1 |
| Built: | 2024-11-26 06:29:35 UTC |
| Source: | CRAN |
Applies the BPS Sampler to a Gaussian target distribution, as detailed in Bouchard-Côté et al, 2017. Assume potential of the form
i.e. a Gaussian with mean vector mu and covariance matrix inv(V)
BPSGaussian(V, mu, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), refresh_rate = 1, unit_velocity = FALSE)BPSGaussian(V, mu, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), refresh_rate = 1, unit_velocity = FALSE)
V |
the inverse covariance matrix (or precision matrix) of the Gaussian target distribution. |
mu |
mean of the Gaussian target distribution |
n_iter |
Number of algorithm iterations; will result in the equivalent amount of skeleton points in Gaussian case because no rejections are needed. |
finalTime |
If provided and nonnegative, run the BPS sampler until a trajectory of continuous time length finalTime is obtained (ignoring the value of |
x0 |
starting point (optional, if not specified taken to be the origin) |
v0 |
starting direction (optional, if not specified taken to be a random vector) |
refresh_rate |
|
unit_velocity |
TRUE indicates velocities uniform on unit sphere, FALSE (default) indicates standard normal velocities |
Returns a list with the following objects:
Times: Vector of switching times
Positions: Matrix whose columns are locations of switches. The number of columns is identical to the length of skeletonTimes. Be aware that the skeleton points themselves are NOT samples from the target distribution.
Velocities: Matrix whose columns are velocities just after switches. The number of columns is identical to the length of skeletonTimes.
V <- matrix(c(3,1,1,3),nrow=2) mu <- c(2,2) x0 <- c(0,0) result <- BPSGaussian(V, mu, n_iter = 100, x0 = x0) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)V <- matrix(c(3,1,1,3),nrow=2) mu <- c(2,2) x0 <- c(0,0) result <- BPSGaussian(V, mu, n_iter = 100, x0 = x0) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)
Applies the Bouncy Particle Sampler to a IID Gaussian distribution
BPSIIDGaussian(variance, dim = 1L, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), refresh_rate = 1, unit_velocity = FALSE)BPSIIDGaussian(variance, dim = 1L, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), refresh_rate = 1, unit_velocity = FALSE)
variance |
scalar indicating variance |
dim |
dimension |
n_iter |
Number of algorithm iterations; will result in the equivalent amount of skeleton points in Gaussian case because no rejections are needed. |
finalTime |
If provided and nonnegative, run the sampler until a trajectory of continuous time length finalTime is obtained (ignoring the value of |
x0 |
starting point (optional, if not specified taken to be the origin) |
v0 |
starting direction (optional, if not specified taken to be a random vector) |
refresh_rate |
|
unit_velocity |
TRUE indicates velocities uniform on unit sphere, FALSE (default) indicates standard normal velocities |
Returns a list with the following objects:
Times: Vector of switching times
Positions: Matrix whose columns are locations of switches. The number of columns is identical to the length of skeletonTimes. Be aware that the skeleton points themselves are NOT samples from the target distribution.
Velocities: Matrix whose columns are velocities just after switches. The number of columns is identical to the length of skeletonTimes.
result <- BPSIIDGaussian(1, 2, 1000) plot(result$Positions[2,], result$Positions[1,],type='l',asp=1)result <- BPSIIDGaussian(1, 2, 1000) plot(result$Positions[2,], result$Positions[1,],type='l',asp=1)
Applies the Zig-Zag Sampler to a Student T distribution (IID or spherically symmetric)
BPSStudentT(dof, dim = 1L, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), sphericallySymmetric = TRUE, refresh_rate = 1, unit_velocity = FALSE)BPSStudentT(dof, dim = 1L, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), sphericallySymmetric = TRUE, refresh_rate = 1, unit_velocity = FALSE)
dof |
scalar indicating degrees of freedom |
dim |
dimension |
n_iter |
Number of algorithm iterations; will result in the equivalent amount of skeleton points in Gaussian case because no rejections are needed. |
finalTime |
If provided and nonnegative, run the sampler until a trajectory of continuous time length finalTime is obtained (ignoring the value of |
x0 |
starting point (optional, if not specified taken to be the origin) |
v0 |
starting direction (optional, if not specified taken to be a random vector) |
sphericallySymmetric |
boolean. If false, sample iid Student T distribution, if true (default) sample spherically summetric Student t dsitribution. |
refresh_rate |
|
unit_velocity |
TRUE indicates velocities uniform on unit sphere, FALSE (default) indicates standard normal velocities |
Returns a list with the following objects:
Times: Vector of switching times
Positions: Matrix whose columns are locations of switches. The number of columns is identical to the length of skeletonTimes. Be aware that the skeleton points themselves are NOT samples from the target distribution.
Velocities: Matrix whose columns are velocities just after switches. The number of columns is identical to the length of skeletonTimes.
dim = 2 dof = 4 result <- BPSStudentT(dof, dim, n_iter=1000,sphericallySymmetric = TRUE) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)dim = 2 dof = 4 result <- BPSStudentT(dof, dim, n_iter=1000,sphericallySymmetric = TRUE) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)
Extract discrete samples from a skeleton
DiscreteSamples(skeletonList, n_samples, coordinate = -1L)DiscreteSamples(skeletonList, n_samples, coordinate = -1L)
skeletonList |
a piecewise deterministic skeleton (consisting of Times, Points and Velocities) returned by a sampler |
n_samples |
number of samples to obtain |
coordinate |
if specified, only obtain samples of the specified coordinate, otherwise obtain samples of all coordinates |
Returns a list containing the extracted samples and the times (on the continuous time scale) at which the samples are extracted
Estimates the covariance matrix of a piecewise deterministic skeleton
EstimateCovarianceMatrix(skeletonList, coordinate = -1L, zeroMeans = FALSE)EstimateCovarianceMatrix(skeletonList, coordinate = -1L, zeroMeans = FALSE)
skeletonList |
a piecewise deterministic skeleton (consisting of Times, Points and Velocities) returned by a sampler |
coordinate |
if specified, only estimate the variance of the specified coordinate, otherwise estimate the covariance matrix of all coordinates |
zeroMeans |
if TRUE do not estimate means but assume a centered distribution |
Returns a list containing the estimated moment
Estimates the effective sample size (ESS) of a piecewise deterministic skeleton
EstimateESS(skeletonList, n_batches = 100L, coordinate = -1L, zeroMeans = FALSE)EstimateESS(skeletonList, n_batches = 100L, coordinate = -1L, zeroMeans = FALSE)
skeletonList |
a piecewise deterministic skeleton (consisting of Times, Points and Velocities) returned by a sampler |
n_batches |
optional argument indicating the number of batches to use in the batch means estimation method |
coordinate |
if specified, only estimate the ESS of the specified coordinate, otherwise estimate the ESS of all coordinates |
zeroMeans |
if TRUE do not estimate means but assume a centered distribution |
Returns a list containing the estimated asymptotic variance, ESS and estimated covariance matrix
Estimates the p-th moment of a piecewise deterministic skeleton
EstimateMoment(skeletonList, p, coordinate = -1L)EstimateMoment(skeletonList, p, coordinate = -1L)
skeletonList |
a piecewise deterministic skeleton (consisting of Times, Points and Velocities) returned by a sampler |
p |
moment to estimate |
coordinate |
if specified, only estimate the ESS of the specified coordinate, otherwise estimate the ESS of all coordinates |
Returns a list containing the estimated moment
Implements various piecewise deterministic Monte Carlo methods, including the Zig-Zag Sampler (Bierkens, Fearnhead, Roberts, 2019, https://arxiv.org/abs/1607.03188) and the Bouncy Particle Sampler (BPS, Bouchard-Côté et al., 2017, https://arxiv.org/abs/1510.02451). Typical usage consists of first creating a "skeleton" consisting of "events", which can be used directly for plotting trajectories. The skeleton may be post-processed to extract information, such as as moment and covariance estimates, discrete samples at fixed time intervals along the trajectory, effective sample size and asymptotic variance.
This package currently consists of the following functions for generating skeletons:
ZigZagLogistic for logistic regression, ZigZagGaussian for multivariate Gaussian, ZigZagIIDGaussian for a IID Gaussian target using Zig-Zag, ZigZagStudentT for spherically symmetric or factorized Student-t distribution, BPSGaussian for multivariate Gaussian using BPS, BPSIIDGaussian for a IID Gaussian target using BPS, BPSStudentT for BPS applied to a spherically symmetric or factorized Student-t distribution.
Furthermore the package contains the following functions for post-processing:
EstimateESS (to estimate asymptotic variance and effective sample size for individual coordinates), EstimateMoment, EstimateCovarianceMatrix and DiscreteSamples.
Joris Bierkens
With thanks to Matt Moores, https://mattstats.wordpress.com/, for his help in getting from C++ code to a CRAN-ready Rcpp based package.
Applies the Zig-Zag Sampler to a Gaussian target distribution, as detailed in Bierkens, Fearnhead, Roberts, The Zig-Zag Process and Super-Efficient Sampling for Bayesian Analysis of Big Data, 2016. Assume potential of the form
i.e. a Gaussian with mean vector mu and covariance matrix inv(V)
ZigZagGaussian(V, mu, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0))ZigZagGaussian(V, mu, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0))
V |
the inverse covariance matrix (or precision matrix) of the Gaussian target distribution. |
mu |
mean of the Gaussian target distribution |
n_iter |
Number of algorithm iterations; will result in the equivalent amount of skeleton points in Gaussian case because no rejections are needed. |
finalTime |
If provided and nonnegative, run the sampler until a trajectory of continuous time length finalTime is obtained (ignoring the value of |
x0 |
starting point (optional, if not specified taken to be the origin) |
v0 |
starting direction (optional, if not specified taken to be +1 in every component) |
Returns a list with the following objects:
Times: Vector of switching times
Positions: Matrix whose columns are locations of switches. The number of columns is identical to the length of skeletonTimes. Be aware that the skeleton points themselves are NOT samples from the target distribution.
Velocities: Matrix whose columns are velocities just after switches. The number of columns is identical to the length of skeletonTimes.
V <- matrix(c(3,1,1,3),nrow=2) mu <- c(2,2) result <- ZigZagGaussian(V, mu, 100) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)V <- matrix(c(3,1,1,3),nrow=2) mu <- c(2,2) result <- ZigZagGaussian(V, mu, 100) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)
Applies the Zig-Zag Sampler to a IID Gaussian distribution
ZigZagIIDGaussian(variance, dim = 1L, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0))ZigZagIIDGaussian(variance, dim = 1L, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0))
variance |
scalar indicating variance |
dim |
dimension |
n_iter |
Number of algorithm iterations; will result in the equivalent amount of skeleton points in Gaussian case because no rejections are needed. |
finalTime |
If provided and nonnegative, run the sampler until a trajectory of continuous time length finalTime is obtained (ignoring the value of |
x0 |
starting point (optional, if not specified taken to be the origin) |
v0 |
starting direction (optional, if not specified taken to be +1 in every component) |
Returns a list with the following objects:
Times: Vector of switching times
Positions: Matrix whose columns are locations of switches. The number of columns is identical to the length of skeletonTimes. Be aware that the skeleton points themselves are NOT samples from the target distribution.
Velocities: Matrix whose columns are velocities just after switches. The number of columns is identical to the length of skeletonTimes.
result <- ZigZagIIDGaussian(1, 2, 1000) plot(result$Positions[2,], result$Positions[1,],type='l',asp=1)result <- ZigZagIIDGaussian(1, 2, 1000) plot(result$Positions[2,], result$Positions[1,],type='l',asp=1)
Applies the Zig-Zag Sampler to logistic regression, as detailed in Bierkens, Fearnhead, Roberts, The Zig-Zag Process and Super-Efficient Sampling for Bayesian Analysis of Big Data, 2019.
ZigZagLogistic(dataX, dataY, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), cv = FALSE)ZigZagLogistic(dataX, dataY, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), cv = FALSE)
dataX |
Design matrix containing observations of the independent variables x. The i-th row represents the i-th observation with components x_i,1, ..., x_i,d. |
dataY |
Vector of length n containing 0, 1-valued observations of the dependent variable y. |
n_iter |
Number of algorithm iterations; will result in the equivalent amount of skeleton points in Gaussian case because no rejections are needed. |
finalTime |
If provided and nonnegative, run the sampler until a trajectory of continuous time length finalTime is obtained (ignoring the value of |
x0 |
starting point (optional, if not specified taken to be the origin) |
v0 |
starting direction (optional, if not specified taken to be +1 in every component) |
cv |
optional boolean to indicate the use of subsampling with control variates |
Returns a list with the following objects:
Times: Vector of switching times
Positions: Matrix whose columns are locations of switches. The number of columns is identical to the length of skeletonTimes. Be aware that the skeleton points themselves are NOT samples from the target distribution.
Velocities: Matrix whose columns are velocities just after switches. The number of columns is identical to the length of skeletonTimes.
require("RZigZag") generate.logistic.data <- function(beta, n.obs) { dim <- length(beta) dataX <- cbind(rep(1.0,n.obs), matrix(rnorm((dim -1) * n.obs), ncol = dim -1)); 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)) } beta <- c(1,2) data <- generate.logistic.data(beta, 1000) result <- ZigZagLogistic(data$dataX, data$dataY, 1000) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)require("RZigZag") generate.logistic.data <- function(beta, n.obs) { dim <- length(beta) dataX <- cbind(rep(1.0,n.obs), matrix(rnorm((dim -1) * n.obs), ncol = dim -1)); 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)) } beta <- c(1,2) data <- generate.logistic.data(beta, 1000) result <- ZigZagLogistic(data$dataX, data$dataY, 1000) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)
Applies the Zig-Zag Sampler to a Student T distribution (IID or spherically symmetric)
ZigZagStudentT(dof, dim = 1L, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), sphericallySymmetric = TRUE)ZigZagStudentT(dof, dim = 1L, n_iter = -1L, finalTime = -1, x0 = numeric(0), v0 = numeric(0), sphericallySymmetric = TRUE)
dof |
scalar indicating degrees of freedom |
dim |
dimension |
n_iter |
Number of algorithm iterations; will result in the equivalent amount of skeleton points in Gaussian case because no rejections are needed. |
finalTime |
If provided and nonnegative, run the sampler until a trajectory of continuous time length finalTime is obtained (ignoring the value of |
x0 |
starting point (optional, if not specified taken to be the origin) |
v0 |
starting direction (optional, if not specified taken to be +1 in every component) |
sphericallySymmetric |
boolean. If false, sample iid Student T distribution, if true (default) sample spherically summetric Student t dsitribution. |
Returns a list with the following objects:
Times: Vector of switching times
Positions: Matrix whose columns are locations of switches. The number of columns is identical to the length of skeletonTimes. Be aware that the skeleton points themselves are NOT samples from the target distribution.
Velocities: Matrix whose columns are velocities just after switches. The number of columns is identical to the length of skeletonTimes.
dim = 2 dof = 4 result <- ZigZagStudentT(dof, dim, n_iter=1000, sphericallySymmetric = TRUE) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)dim = 2 dof = 4 result <- ZigZagStudentT(dof, dim, n_iter=1000, sphericallySymmetric = TRUE) plot(result$Positions[1,], result$Positions[2,],type='l',asp=1)