ZigZagLogistic {RZigZag}

ZigZagLogistic

Description

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.

Usage

ZigZagLogistic(dataX, dataY, n_iter = -1L, finalTime = -1,
  x0 = numeric(0), v0 = numeric(0), cv = FALSE)

Arguments

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

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

Value

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.

Examples

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)