# "Shake and Bake" sampler 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. ```r 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 - `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`. `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. ## Returns 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`. - **x0**: The generated seed point. See `createSeedPoint`. - **i0**: 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). ## Author(s) Gert van Valkenhoef ## Note "Shake and Bake" is a Markov Chain Monte Carlo (MCMC) method, so generated samples form a correlated time series. ## See Also `harConstraints` `sab` ## Examples ```r # 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) ```

shakeandbake function