hitandrun: "Hit and Run" sampler

In hitandrun: "Hit and Run" and "Shake and Bake" for Sampling Uniformly from Convex Shapes

"Hit and Run" sampler

Description

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

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

Usage

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

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

har.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
x0	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)

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 ≤q b (depending on dir). The lengths of rhs and dir must match the number of rows of constr.

hitandrun 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. Finally, har is used to generate the samples.

Value

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

For har.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.
x0	The generated seed point. See createSeedPoint.
thin	The thinning factor to be used.

For har.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 has been modified).

Note

"Hit and Run" is a Markov Chain Monte Carlo (MCMC) method, so generated samples form a correlated time series. To get a uniform sample, you need O*(n^3) samples, where n is the dimension of the sampling space.

Author(s)

Gert van Valkenhoef

See Also

harConstraints har

Examples

# 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 <- hitandrun(constr, n.samples=1000)
stopifnot(dim(samples) == c(1000, 3))
stopifnot(all.equal(apply(samples, 1, sum), rep(1, 1000)))
stopifnot(all.equal(samples[,1]/samples[,2], rep(2, 1000)))

# 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 <- har.init(constr)
result <- har.run(state, n.samples=1000)
samples <- result$samples
stopifnot(all(samples >= 0 & samples <= 1))
# Continue sampling from the same chain:
result <- har.run(result$state, n.samples=1000)
samples <- rbind(samples, result$samples)

hitandrun documentation built on May 28, 2022, 1:09 a.m.