Markov Chain Monte Carlo Sampling#

Markov Chains#

Proposal Distributions#

hopsy is not particulary designed for development of polytope samplers, but instead aims at making the proposals available in hops easily usable to solve the practitioner’s problems at hand. As the proposal - together with the Metropolis filter and the models likelihood function - forms the core of any Markov chain Monte Carlo algorithm, we believe that for the sake of performance proposals should be optimized and implemented in C++ and hops itself.

Nevertheless, it might be desirable to do some rapid prototyping to test promising ideas. For this case, it is possible to implement Python proposals and to instruct hopsy (and ultimately hops) to use those.

Similarily to the custom Python models, this works by wrapping the Python proposal in the hopsy.PyProposal class. hopsy.PyProposal implements all functions necessary for the proposal class to be usable in hops by delegating the function calls to the stored Python proposal. Thus, it is obviously needed, that the Python proposal implements the corresponding functions.

The required functions are

Note

THIS SECTION IS DEPRECATED AND NEEDS TO BE REWRITTEN.

  • propose() -> None generates a new proposal and store it internally

  • accept_proposal() -> None sets the current state, which also has to be stored internally, to the current proposal

  • calculate_log_acceptance_probability() -> float compute the logarithm of \(\frac{P(x|y)}{P(y|x)}\), if \(x\) is the current state and \(y\) the proposal. For symmetric proposals this is e.g. 0. For infeasible proposals, which e.g. lie outside the polytope, it is \(-\infty\).

  • get_state() -> numpy.ndarray[shape[d,1]] returns the current state.

  • set_state(new_state: numpy.ndarray[shape[d,1]]) -> None sets the current state.

  • get_proposal() -> numpy.ndarray[shape[d,1]] returns the current proposal.

  • (optional) get_stepsize() -> float sets the stepsize, if it is available.

  • (optional) set_stepsize(new_stepsize: float) gets the stepsize, if it is available.

  • (optional) get_name() -> str gets the algorithms name, if it is available.

Example code#

Reference#

Proposals#

hopsy.AdaptiveMetropolisProposal

Adaptive Metropolis Algorithm for convex polytopes, see https://doi.org/10.2307/3318737 for unconstrained version.

hopsy.BallWalkProposal

This proposal is implemented for academic purpose and not recommended for production use, see https://doi.org/10.1007/s12532-015-0097-z for description.

hopsy.CSmMALAProposal

Do not use when fisher information is expensive for your model.

hopsy.DikinWalkProposal

See https://doi.org/10.1016/j.orl.2016.07.005.

hopsy.GaussianCoordinateHitAndRunProposal

Mix of https://doi.org/10.1093/bioinformatics/btx052 and https://deepblue.lib.umich.edu/bitstream/handle/2027.42/3513/bal7884.0001.001.pdf.

hopsy.GaussianHitAndRunProposal

Mix of https://doi.org/10.1093/bioinformatics/btx052 and https://deepblue.lib.umich.edu/bitstream/handle/2027.42/3513/bal7884.0001.001.pdf.

hopsy.GaussianProposal

Random walk with gaussian steps.

hopsy.UniformCoordinateHitAndRunProposal

See https://doi.org/10.1093/bioinformatics/btx052.

hopsy.UniformHitAndRunProposal

See: https://dl.acm.org/doi/pdf/10.1145/256562.256619

Markov Chain#

hopsy.MarkovChain(problem[, proposal, ...])

Given a hopsy.Problem a MarkovChain object can be constructed.

hopsy.sample(markov_chains, rngs, n_samples)

Draw n_samples from every passed chain in markov_chains using the respective random number generator from rngs.

Random#

hopsy.RandomNumberGenerator

hopsy.Normal

hopsy.Uniform