shakeandbake: "Shake and Bake" sampler

Description

The "Shake and Bake" method generates a Markov Chain whose stable state converges on the uniform distribution over a the boundary of a convex polytope defined by a set of linear inequality constraints. shakeandbake further uses the Moore-Penrose pseudo-inverse to eliminate an arbitrary set of linear equality constraints before applying the "Shake and Bake" sampler.

sab.init and sab.run together provide a re-entrant version of shakeandbake so that the Markov chain can be continued if convergence is not satisfactory.

Usage

shakeandbake(constr, n.samples=1E4,
    thin.fn = function(n) { ceiling(log(n + 1)/4 * n^3) }, thin = NULL,
    x0.randomize=FALSE, x0.method="slacklp", x0 = NULL, eliminate = TRUE)
sab.init(constr,
    thin.fn = function(n) { ceiling(log(n + 1)/4 * n^3) }, thin = NULL,
    x0.randomize=FALSE, x0.method="slacklp", x0 = NULL, eliminate = TRUE)
sab.run(state, n.samples)

Arguments

constr

Linear constraints that define the sampling space (see details)

n.samples

The desired number of samples to return. The sampler is run for n.samples * thin iterations

thin.fn

Function that specifies a thinning factor depending on the dimension of the sampling space after equality constraints have been eliminated. Will only be invoked if thin is NULL

thin

The thinning factor

Seed point for the Markov Chain. The seed point is specified in the original space, and transformed to the sampling space automatically.

x0.method

Method to generate the seed point if x0 is unspecified, see createSeedPoint

x0.randomize

Whether to generate a random seed point if x0 is unspecified

eliminate

Whether to eliminate redundant constraints before constructing the transformation to the sampling space and (optionally) calculating the seed point.

state

A state object, as generated by har.init (see value)

Value

For shakeandbake, a matrix containing the generated samples as rows.

For sab.init, a state object, containing:

basis

The basis for the sampling space. See solution.basis.

transform

The sampling space transformation. See createTransform.

constr

The linear inequality constraints translated to the sampling space. See transformConstraints.

The generated seed point. See createSeedPoint.

The index of the closest face. See findFace.

thin

The thinning factor to be used.

For sab.run, a list containing:

samples

A matrix containing the generated samples as rows.

state

A state object that can be used to continue sampling from the Markov chain (i.e. x0 and i0 have been modified).

Details

The constraints are given as a list with the elements constr, dir and rhs. dir is a vector with values '=' or '<='. constr is a matrix and rhs a vector, which encode the standard linear programming constraint froms \(Ax = b\) and \(Ax \leq b\) (depending on dir). The lengths of rhs and dir must match the number of rows of constr.

shakeandbake applies solution.basis to generate a basis of the (translated) solution space of the linear constraints (if any). An affine transformation is generated using createTransform and applied to the constraints. Then, a seed point satisfying the inequality constraints is generated using createSeedPoint. The closest face to this point is found using findFace. Finally, sab is used to generate the samples.

Examples

Run this code

# NOT RUN {
# Sample from the 3-simplex with the additional constraint that w_1/w_2 = 2
# Three inequality constraints, two equality constraints
constr <- mergeConstraints(simplexConstraints(3), exactRatioConstraint(3, 1, 2, 2))
samples <- shakeandbake(constr, n.samples=1000)
stopifnot(dim(samples) == c(1000, 3))
stopifnot(all.equal(apply(samples, 1, sum), rep(1, 1000)))

sel <- samples[,3] > 0.5 # detect which side we're on
stopifnot(all.equal(samples[sel,], matrix(rep(c(0,0,1), each=sum(sel)), ncol=3)))
stopifnot(all.equal(samples[!sel,], matrix(rep(c(2/3,1/3,0), each=sum(sel)), ncol=3)))

# Sample from the unit rectangle (no equality constraints)
constr <- list(
  constr = rbind(c(1,0), c(0,1), c(-1,0), c(0,-1)),
  dir=rep('<=', 4),
  rhs=c(1, 1, 0, 0))
state <- sab.init(constr)
result <- sab.run(state, n.samples=1000)
faces <- result$faces
samples <- result$samples
stopifnot(all(samples >= -1e-15 & samples <= 1 + 1e-15))

stopifnot(all.equal(samples[faces==1,1], rep(1, sum(faces==1))))
stopifnot(all.equal(samples[faces==2,2], rep(1, sum(faces==2))))
stopifnot(all.equal(samples[faces==3,1], rep(0, sum(faces==3))))
stopifnot(all.equal(samples[faces==4,2], rep(0, sum(faces==4))))

# Continue sampling from the same chain:
result <- sab.run(result$state, n.samples=1000)
samples <- rbind(samples, result$samples)
# }