hitandrun {hitandrun}

"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 \leq 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)