mici.samplers module#

Monte Carlo sampler classes for peforming inference.

class mici.samplers.DynamicMultinomialHMC(system, integrator, rng, *, max_tree_depth=10, max_delta_h=1000, termination_criterion=<function riemannian_no_u_turn_criterion>, do_extra_subtree_checks=True, momentum_transition=None)[source]#

Bases: HamiltonianMonteCarlo

Dynamic integration time HMC method with multinomial sampling from trajectory.

In each transition a binary tree of states is recursively computed by integrating randomly forward and backward in time by a number of steps equal to the previous tree size (Hoffman and Gelman, 2014; Betancourt, 2017) until a termination criteria on the tree leaves is met. The next chain state is chosen from the candidate states using a progressive multinomial sampling scheme (Betancourt, 2017) based on the relative probability densities of the different candidate states, with the resampling biased towards states further from the current state.

When used with the default settings of mici.transitions.riemannian_no_u_turn_criterion termination criterion and extra subtree checks enabled (do_extra_subtree_checks == True), this sampler is equivalent to the default ‘NUTS’ MCMC algorithm (minus adaptation) used in Stan as of version v2.23.

References

Hoffman, M.D. and Gelman, A., 2014. The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1), pp.1593-1623.
Betancourt, M., 2017. A conceptual introduction to Hamiltonian Monte Carlo. arXiv preprint arXiv:1701.02434.

Parameters:

system (System) – Hamiltonian system to be simulated.
rng (Generator) – Numpy random number generator.
integrator (Integrator) – Symplectic integrator to use to simulate dynamics in integration transition.
max_tree_depth (int) – Maximum depth to expand trajectory binary tree to in integrator transition. The maximum number of integrator steps corresponds to 2**max_tree_depth.
max_delta_h (float) – Maximum change to tolerate in the Hamiltonian function over a trajectory in integrator transition before signalling a divergence.
termination_criterion (TerminationCriterion) – Function computing criterion to use to determine when to terminate trajectory tree expansion. The function should take a Hamiltonian system as its first argument, a pair of states corresponding to the two edge nodes in the trajectory (sub-)tree being checked and an array containing the sum of the momentums over the trajectory (sub)-tree. Defaults to mici.transitions.riemannian_no_u_turn_criterion.
do_extra_subtree_checks (bool) – Whether to perform additional termination criterion checks on overlapping subtrees of the current tree to improve robustness in systems with dynamics which are well approximated by independent system of simple harmonic oscillators. In such systems (corresponding to e.g. a standard normal target distribution and identity metric matrix representation) at certain step sizes a ‘resonant’ behaviour is seen by which the termination criterion fails to detect that the trajectory has expanded past a half-period i.e. has ‘U-turned’ resulting in trajectories continuing to expand, potentially up until the max_tree_depth limit is hit. For more details see this Stan Discourse discussion. If do_extra_subtree_checks is set to True additional termination criterion checks are performed on overlapping subtrees which help to reduce this resonant behaviour at the cost of more conservative trajectory termination in some correlated models and some overhead from additional checks.
momentum_transition (Optional[MomentumTransition]) – Markov transition kernel which leaves the conditional distribution on the momentum under the canonical distribution invariant, updating only the momentum component of the chain state. If set to None the momentum transition operator mici.transitions.IndependentMomentumTransition will be used, which independently samples the momentum from its conditional distribution.

property max_delta_h: float#: Change in Hamiltonian over trajectory to trigger divergence.

property max_tree_depth: int#: Maximum depth to expand trajectory binary tree to.

property rng: Generator#: NumPy random number generator object.

sample_chains(n_warm_up_iter, n_main_iter, init_states, **kwargs)#

Sample Markov chains from given initial states with adaptive warm up.

One or more Markov chains are sampled, with each chain iteration consisting of a momentum transition followed by an integration transition. The chains are split into multiple stages with zero or more adaptive warm up stages followed by the main non-adaptive sampling stage. During the adaptive stage(s) parameters of the integration transition such as the integrator step size are adaptively tuned based on the chain state and/or transition statistics.

The default settings use a single (fast) DualAveragingStepSizeAdapter adapter instance which adapts the integrator step-size using a dual-averaging algorithm in a single adaptive stage.

The chains (including both adaptive and non-adaptive stages) may be run in parallel across multiple independent processes or sequentially. In all cases all chains use independent random draws.

Parameters:

n_warm_up_iter (int) – Number of adaptive warm up iterations per chain. Depending on the mici.stagers.Stager instance specified by the stager argument the warm up iterations may be split between one or more adaptive stages. If zero, only a single non-adaptive stage is used.
n_main_iter (int) – Number of iterations (samples to draw) per chain during main (non-adaptive) sampling stage.
init_states (Iterable[Union[ChainState, NDArray, dict]]) – Initial chain states. Each state can be either an array specifying the state position component or a mici.states.ChainState instance. If an array is passed or the mom attribute of the state is not set, a momentum component will be independently sampled from its conditional distribution. One chain will be run for each state in the iterable.

Keyword Arguments:

trace_funcs – Sequence of functions which compute the variables to be recorded at each chain iteration (during only the main non-adaptive sampling stage if trace_warm_up is False), with each trace function passed the current state and returning a dictionary of scalar or array values corresponding to the variable(s) to be stored. The keys in the returned dictionaries are used to index the trace arrays in the returned traces dictionary. If a key appears in multiple dictionaries only the the value corresponding to the last trace function to return that key will be stored. Default is to use a single function which records the position component of the state under the key "pos" and the Hamiltonian at the state under the key "hamiltonian".
adapters – Sequence of mici.adapters.Adapter instances to use to adaptatively set parameters of the integration transition such as the step size during the adaptive stages of the chains. Note that the adapter updates are applied in the order the adapters appear in the sequence and so if multiple adapters change the same parameter(s) the order will matter. If None or an empty sequence no adapters are used. Default is to use a single instance of mici.adapters.DualAveragingStepSizeAdapter with its default parameters.
stager – Chain iteration stager object which controls the split of the chain iterations into the adaptive warm up and non-adaptive main stages. If set to None (the default) and all adapters specified by the adapters argument are of the fast type (i.e. their is_fast attribute is True) then a mici.stagers.WarmUpStager instance will be used corresponding to using a single adaptive warm up stage will all adapters active. If set to None and the adapters specified by the adapters argument are not all of the fast type, then a mici.stagers.WindowedWarmUpStager (with its default arguments) will be used, corresponding to using multiple adaptive warm up stages with only the fast-type adapters active in some - see documentation of mici.stagers.WarmUpStager for details.
n_process – Number of parallel processes to run chains over. If n_process=1 then chains will be run sequentially otherwise a multiprocessing.Pool object will be used to dynamically assign the chains across multiple processes. If set to None then the number of processes will be set to the output of os.cpu_count(). Default is n_process=1.
trace_warm_up – Whether to record chain traces and statistics during warm-up stage iterations (True) or only record traces and statistics in the iterations of the final non-adaptive stage (False, the default).
max_threads_per_process – If threadpoolctl is available this argument may be used to limit the maximum number of threads that can be used in thread pools used in libraries supported by threadpoolctl, which include BLAS and OpenMP implementations. This argument will only have an effect if n_process > 1 such that chains are being run on multiple processes and only if threadpoolctl is installed in the current Python environment. If set to None (the default) no limits are set.
force_memmap – Whether to force arrays used to store chain data to be memory-mapped to files on disk to avoid excessive system memory usage for long chains and/or large chain states. The chain data is written to .npy files in the directory specified by memmap_path (or a temporary directory if memmap_path is None). Chain data is always memory mapped when sampling chains in parallel on multiple processes.
memmap_path – Path to directory to write memory-mapped chain data to. If None (the default) and memory-mapping is enabled then a temporary directory will be created and the chain data written to files there, with the created files being deleted in this case once the last reference to them is closed.
monitor_stats – Sequence of string keys of (integration) transition statistics to monitor mean of over samples computed so far during sampling by printing as postfix to progress bar. Default is to print only the mean accept_stat (acceptance statistic).
display_progress – Whether to display a progress bar to track the completed chain sampling iterations. Default value is True, i.e. to display progress bar.
progress_bar_class – Class or factory function for progress bar to use to show chain progress if enabled (display_progress=True). Defaults to mici.progressbars.SequenceProgressBar.

Returns:

Named tuple (final_states, traces, statistics) corresponding to states of chains after final iteration, dictionary of chain trace arrays and dictionary of chain statistics dictionaries.

Return type:

HMCSampleChainsOutputs

property system: System#: Hamiltonian system corresponding to joint distribution on augmented space.

property transitions: dict[str, Transition]#: Dictionary of transition kernels sampled from in each chain iteration.

class mici.samplers.DynamicSliceHMC(system, integrator, rng, *, max_tree_depth=10, max_delta_h=1000.0, termination_criterion=<function euclidean_no_u_turn_criterion>, do_extra_subtree_checks=False, momentum_transition=None)[source]#

Bases: HamiltonianMonteCarlo

Dynamic integration time HMC method with slice sampling from trajectory.

In each transition a binary tree of states is recursively computed by integrating randomly forward and backward in time by a number of steps equal to the previous tree size (Hoffman and Gelman, 2014) until a termination criteria on the tree leaves is met. The next chain state is chosen from the candidate states using a progressive slice sampling scheme (Hoffman and Gelman, 2014) based on the relative probability densities of the different candidate states, with the sampling biased towards states further from the current state.

When used with the default setting of euclidean_no_u_turn_criterion termination criterion and extra subtree checks disabled, this sampler is equivalent to ‘Algorithm 3: Efficient No-U-Turn Sampler’ in Hoffman and Gelman (2014), i.e. the ‘classic NUTS’ algorithm.

References

Hoffman, M.D. and Gelman, A., 2014. The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1), pp.1593-1623.

Parameters:

system (System) – Hamiltonian system to be simulated.
rng (Generator) – Numpy random number generator.
integrator (Integrator) – Symplectic integrator to use to simulate dynamics in integration transition.
max_tree_depth (int) – Maximum depth to expand trajectory binary tree to in integrator transition. The maximum number of integrator steps corresponds to 2**max_tree_depth.
max_delta_h (float) – Maximum change to tolerate in the Hamiltonian function over a trajectory in integrator transition before signalling a divergence.
termination_criterion (TerminationCriterion) – Function computing criterion to use to determine when to terminate trajectory tree expansion. The function should take a Hamiltonian system as its first argument, a pair of states corresponding to the two edge nodes in the trajectory (sub-)tree being checked and an array containing the sum of the momentums over the trajectory (sub)-tree. Defaults to mici.transitions.euclidean_no_u_turn_criterion.
do_extra_subtree_checks (bool) –
Whether to perform additional termination criterion checks on overlapping subtrees of the current tree to improve robustness in systems with dynamics which are well approximated by independent system of simple harmonic oscillators. In such systems (corresponding to e.g. a standard normal target distribution and identity metric matrix representation) at certain step sizes a ‘resonant’ behaviour is seen by which the termination criterion fails to detect that the trajectory has expanded past a half-period i.e. has ‘U-turned’ resulting in trajectories continuing to expand, potentially up until the max_tree_depth limit is hit. For more details see this Stan Discourse discussion. If do_extra_subtree_checks is set to True additional termination criterion checks are performed on overlapping subtrees which help to reduce this resonant behaviour at the cost of more conservative trajectory termination in some correlated models and some overhead from additional checks.
momentum_transition (Optional[MomentumTransition]) – Markov transition kernel which leaves the conditional distribution on the momentum under the canonical distribution invariant, updating only the momentum component of the chain state. If set to None the momentum transition operator mici.transitions.IndependentMomentumTransition will be used, which independently samples the momentum from its conditional distribution.

property max_delta_h: float#: Change in Hamiltonian over trajectory to trigger divergence.

property max_tree_depth: int#: Maximum depth to expand trajectory binary tree to.

property rng: Generator#: NumPy random number generator object.

sample_chains(n_warm_up_iter, n_main_iter, init_states, **kwargs)#

Sample Markov chains from given initial states with adaptive warm up.

The default settings use a single (fast) DualAveragingStepSizeAdapter adapter instance which adapts the integrator step-size using a dual-averaging algorithm in a single adaptive stage.

The chains (including both adaptive and non-adaptive stages) may be run in parallel across multiple independent processes or sequentially. In all cases all chains use independent random draws.

Parameters:

n_warm_up_iter (int) – Number of adaptive warm up iterations per chain. Depending on the mici.stagers.Stager instance specified by the stager argument the warm up iterations may be split between one or more adaptive stages. If zero, only a single non-adaptive stage is used.
n_main_iter (int) – Number of iterations (samples to draw) per chain during main (non-adaptive) sampling stage.
init_states (Iterable[Union[ChainState, NDArray, dict]]) – Initial chain states. Each state can be either an array specifying the state position component or a mici.states.ChainState instance. If an array is passed or the mom attribute of the state is not set, a momentum component will be independently sampled from its conditional distribution. One chain will be run for each state in the iterable.

Keyword Arguments:

trace_funcs – Sequence of functions which compute the variables to be recorded at each chain iteration (during only the main non-adaptive sampling stage if trace_warm_up is False), with each trace function passed the current state and returning a dictionary of scalar or array values corresponding to the variable(s) to be stored. The keys in the returned dictionaries are used to index the trace arrays in the returned traces dictionary. If a key appears in multiple dictionaries only the the value corresponding to the last trace function to return that key will be stored. Default is to use a single function which records the position component of the state under the key "pos" and the Hamiltonian at the state under the key "hamiltonian".
adapters – Sequence of mici.adapters.Adapter instances to use to adaptatively set parameters of the integration transition such as the step size during the adaptive stages of the chains. Note that the adapter updates are applied in the order the adapters appear in the sequence and so if multiple adapters change the same parameter(s) the order will matter. If None or an empty sequence no adapters are used. Default is to use a single instance of mici.adapters.DualAveragingStepSizeAdapter with its default parameters.
stager – Chain iteration stager object which controls the split of the chain iterations into the adaptive warm up and non-adaptive main stages. If set to None (the default) and all adapters specified by the adapters argument are of the fast type (i.e. their is_fast attribute is True) then a mici.stagers.WarmUpStager instance will be used corresponding to using a single adaptive warm up stage will all adapters active. If set to None and the adapters specified by the adapters argument are not all of the fast type, then a mici.stagers.WindowedWarmUpStager (with its default arguments) will be used, corresponding to using multiple adaptive warm up stages with only the fast-type adapters active in some - see documentation of mici.stagers.WarmUpStager for details.
n_process – Number of parallel processes to run chains over. If n_process=1 then chains will be run sequentially otherwise a multiprocessing.Pool object will be used to dynamically assign the chains across multiple processes. If set to None then the number of processes will be set to the output of os.cpu_count(). Default is n_process=1.
trace_warm_up – Whether to record chain traces and statistics during warm-up stage iterations (True) or only record traces and statistics in the iterations of the final non-adaptive stage (False, the default).
max_threads_per_process – If threadpoolctl is available this argument may be used to limit the maximum number of threads that can be used in thread pools used in libraries supported by threadpoolctl, which include BLAS and OpenMP implementations. This argument will only have an effect if n_process > 1 such that chains are being run on multiple processes and only if threadpoolctl is installed in the current Python environment. If set to None (the default) no limits are set.
force_memmap – Whether to force arrays used to store chain data to be memory-mapped to files on disk to avoid excessive system memory usage for long chains and/or large chain states. The chain data is written to .npy files in the directory specified by memmap_path (or a temporary directory if memmap_path is None). Chain data is always memory mapped when sampling chains in parallel on multiple processes.
memmap_path – Path to directory to write memory-mapped chain data to. If None (the default) and memory-mapping is enabled then a temporary directory will be created and the chain data written to files there, with the created files being deleted in this case once the last reference to them is closed.
monitor_stats – Sequence of string keys of (integration) transition statistics to monitor mean of over samples computed so far during sampling by printing as postfix to progress bar. Default is to print only the mean accept_stat (acceptance statistic).
display_progress – Whether to display a progress bar to track the completed chain sampling iterations. Default value is True, i.e. to display progress bar.
progress_bar_class – Class or factory function for progress bar to use to show chain progress if enabled (display_progress=True). Defaults to mici.progressbars.SequenceProgressBar.

Returns:

Named tuple (final_states, traces, statistics) corresponding to states of chains after final iteration, dictionary of chain trace arrays and dictionary of chain statistics dictionaries.

Return type:

HMCSampleChainsOutputs

property system: System#: Hamiltonian system corresponding to joint distribution on augmented space.

property transitions: dict[str, Transition]#: Dictionary of transition kernels sampled from in each chain iteration.

class mici.samplers.HMCSampleChainsOutputs(final_states, traces, statistics)[source]#

Bases: NamedTuple

Outputs returned by HamiltonianMonteCarlo.sample_chains() call.

Parameters:

final_states (list[ChainState]) – States of chains after final iteration. May be used to resume sampling a chain by passing as the initial states to a new HamiltonianMonteCarlo.sample_chains() call.
traces (dict[str, list[NDArray]]) – Dictionary of chain trace arrays. Values in dictionary are list of arrays of variables outputted by trace functions in trace_funcs with each array in the list corresponding to a single chain and the leading dimension of each array corresponding to the iteration (draw) index, within the main non-adaptive sampling stage if trace_warm_up=False and across both warm-up and main sampling stages otherwise. The key for each value is the corresponding key in the dictionary returned by the trace function which computed the traced value.
statistics (dict[str, list[NDArray]]) – Dictionary of chain transition statistic dictionaries. Values in dictionary are lists of arrays of chain statistic values with each array in the list corresponding to a single chain and the leading dimension of each array corresponding to the iteration (draw) index, within the main non-adaptive sampling stage if trace_warm_up=False and across both warm-up and main sampling stages otherwise. The key for each value is a string description of the corresponding integration transition statistic.

class mici.samplers.HamiltonianMonteCarlo(system, rng, integration_transition, momentum_transition=None)[source]#

Bases: MarkovChainMonteCarloMethod

Wrapper class for Hamiltonian Monte Carlo (HMC) methods.

Here HMC is defined as a Markov chain Monte Carlo method which augments the original target variable (henceforth position variable) with a momentum variable with a user specified conditional distribution given the position variable. In each chain iteration two Markov transitions leaving the resulting joint distribution on position and momentum variables invariant are applied - the momentum variables are updated in a transition which leaves their conditional distribution invariant (momentum transition) and then a trajectory in the joint space is generated by numerically integrating a Hamiltonian dynamic with an appropriate symplectic integrator which is exactly reversible, volume preserving and approximately conserves the joint probability density of the target-momentum state pair; one state from the resulting trajectory is then selected as the next joint chain state using an appropriate sampling scheme such that the joint distribution is left exactly invariant (integration transition).

There are various options available for both the momentum transition and integration transition, with by default the momentum transition set to be independent resampling of the momentum variables from their conditional distribution.

References

Duane, S., Kennedy, A.D., Pendleton, B.J. and Roweth, D., 1987. Hybrid Monte Carlo. Physics letters B, 195(2), pp.216-222.
Neal, R.M., 2011. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2(11), p.2.

Parameters:

system (System) – Hamiltonian system to be simulated, corresponding to joint distribution on augmented space.
rng (Generator) – Numpy random number generator.
integration_transition (IntegrationTransition) – Markov transition kernel which leaves joint distribution invariant and jointly updates the position and momentum components of the chain state by integrating the Hamiltonian dynamics of the system to propose new values for the state.
momentum_transition (Optional[MomentumTransition]) – Markov transition kernel which leaves the conditional distribution on the momentum under the join distribution invariant, updating only the momentum component of the chain state. If set to None the momentum transition operator mici.transitions.IndependentMomentumTransition will be used, which independently samples the momentum from its conditional distribution.

property rng: Generator#: NumPy random number generator object.

sample_chains(n_warm_up_iter, n_main_iter, init_states, **kwargs)[source]#

Sample Markov chains from given initial states with adaptive warm up.

The default settings use a single (fast) DualAveragingStepSizeAdapter adapter instance which adapts the integrator step-size using a dual-averaging algorithm in a single adaptive stage.

The chains (including both adaptive and non-adaptive stages) may be run in parallel across multiple independent processes or sequentially. In all cases all chains use independent random draws.

Parameters:

n_warm_up_iter (int) – Number of adaptive warm up iterations per chain. Depending on the mici.stagers.Stager instance specified by the stager argument the warm up iterations may be split between one or more adaptive stages. If zero, only a single non-adaptive stage is used.
n_main_iter (int) – Number of iterations (samples to draw) per chain during main (non-adaptive) sampling stage.
init_states (Iterable[Union[ChainState, NDArray, dict]]) – Initial chain states. Each state can be either an array specifying the state position component or a mici.states.ChainState instance. If an array is passed or the mom attribute of the state is not set, a momentum component will be independently sampled from its conditional distribution. One chain will be run for each state in the iterable.

Keyword Arguments:

trace_funcs – Sequence of functions which compute the variables to be recorded at each chain iteration (during only the main non-adaptive sampling stage if trace_warm_up is False), with each trace function passed the current state and returning a dictionary of scalar or array values corresponding to the variable(s) to be stored. The keys in the returned dictionaries are used to index the trace arrays in the returned traces dictionary. If a key appears in multiple dictionaries only the the value corresponding to the last trace function to return that key will be stored. Default is to use a single function which records the position component of the state under the key "pos" and the Hamiltonian at the state under the key "hamiltonian".
adapters – Sequence of mici.adapters.Adapter instances to use to adaptatively set parameters of the integration transition such as the step size during the adaptive stages of the chains. Note that the adapter updates are applied in the order the adapters appear in the sequence and so if multiple adapters change the same parameter(s) the order will matter. If None or an empty sequence no adapters are used. Default is to use a single instance of mici.adapters.DualAveragingStepSizeAdapter with its default parameters.
stager – Chain iteration stager object which controls the split of the chain iterations into the adaptive warm up and non-adaptive main stages. If set to None (the default) and all adapters specified by the adapters argument are of the fast type (i.e. their is_fast attribute is True) then a mici.stagers.WarmUpStager instance will be used corresponding to using a single adaptive warm up stage will all adapters active. If set to None and the adapters specified by the adapters argument are not all of the fast type, then a mici.stagers.WindowedWarmUpStager (with its default arguments) will be used, corresponding to using multiple adaptive warm up stages with only the fast-type adapters active in some - see documentation of mici.stagers.WarmUpStager for details.
n_process – Number of parallel processes to run chains over. If n_process=1 then chains will be run sequentially otherwise a multiprocessing.Pool object will be used to dynamically assign the chains across multiple processes. If set to None then the number of processes will be set to the output of os.cpu_count(). Default is n_process=1.
trace_warm_up – Whether to record chain traces and statistics during warm-up stage iterations (True) or only record traces and statistics in the iterations of the final non-adaptive stage (False, the default).
max_threads_per_process – If threadpoolctl is available this argument may be used to limit the maximum number of threads that can be used in thread pools used in libraries supported by threadpoolctl, which include BLAS and OpenMP implementations. This argument will only have an effect if n_process > 1 such that chains are being run on multiple processes and only if threadpoolctl is installed in the current Python environment. If set to None (the default) no limits are set.
force_memmap – Whether to force arrays used to store chain data to be memory-mapped to files on disk to avoid excessive system memory usage for long chains and/or large chain states. The chain data is written to .npy files in the directory specified by memmap_path (or a temporary directory if memmap_path is None). Chain data is always memory mapped when sampling chains in parallel on multiple processes.
memmap_path – Path to directory to write memory-mapped chain data to. If None (the default) and memory-mapping is enabled then a temporary directory will be created and the chain data written to files there, with the created files being deleted in this case once the last reference to them is closed.
monitor_stats – Sequence of string keys of (integration) transition statistics to monitor mean of over samples computed so far during sampling by printing as postfix to progress bar. Default is to print only the mean accept_stat (acceptance statistic).
display_progress – Whether to display a progress bar to track the completed chain sampling iterations. Default value is True, i.e. to display progress bar.
progress_bar_class – Class or factory function for progress bar to use to show chain progress if enabled (display_progress=True). Defaults to mici.progressbars.SequenceProgressBar.

Returns:

Named tuple (final_states, traces, statistics) corresponding to states of chains after final iteration, dictionary of chain trace arrays and dictionary of chain statistics dictionaries.

Return type:

HMCSampleChainsOutputs

property system: System#: Hamiltonian system corresponding to joint distribution on augmented space.

property transitions: dict[str, Transition]#: Dictionary of transition kernels sampled from in each chain iteration.

class mici.samplers.MCMCSampleChainsOutputs(final_states, traces, statistics)[source]#

Bases: NamedTuple

Outputs returned by MarkovChainMonteCarloMethod.sample_chains() call.

Parameters:

final_states (list[ChainState]) – States of chains after final iteration. May be used to resume sampling a chain by passing as the initial states to a new MarkovChainMonteCarloMethod.sample_chains() call.
traces (dict[str, list[NDArray]]) – Dictionary of chain trace arrays. Values in dictionary are list of arrays of variables outputted by trace functions in trace_funcs with each array in the list corresponding to a single chain and the leading dimension of each array corresponding to the iteration (draw) index, within the main non-adaptive sampling stage if trace_warm_up=False and across both warm-up and main sampling stages otherwise. The key for each value is the corresponding key in the dictionary returned by the trace function which computed the traced value.
statistics (dict[str, dict[str, list[NDArray]]]) – Dictionary of chain transition statistic dictionaries. Values in outer dictionary are dictionaries of statistics for each chain transition, keyed by the string key for the transition. The values in each inner transition dictionary are lists of arrays of chain statistic values with each array in the list corresponding to a single chain and the leading dimension of each array corresponding to the iteration (draw) index, within the main non-adaptive sampling stage if trace_warm_up=False and across both warm-up and main sampling stages otherwise. The key for each value is a string description of the corresponding transition statistic.

class mici.samplers.MarkovChainMonteCarloMethod(rng, transitions)[source]#

Bases: object

Generic Markov chain Monte Carlo (MCMC) sampler.

Generates a Markov chain from some initial state by iteratively applying a sequence of Markov transition operators.

Parameters:

rng (Generator) – Numpy random number generator.
transitions (dict[str, Transition]) – Ordered dictionary of Markov transitions kernels to sequentially sample from on each chain iteration.

property rng: Generator#: NumPy random number generator object.

sample_chains(n_warm_up_iter, n_main_iter, init_states, *, trace_funcs=None, adapters=None, stager=None, n_process=1, trace_warm_up=False, max_threads_per_process=None, force_memmap=False, memmap_path=None, monitor_stats=None, display_progress=True, progress_bar_class=None)[source]#

Sample Markov chains from initial states with optional adaptive warm up.

One or more Markov chains are sampled, with each chain iteration consisting of one or more Markov transitions. The chains are split into multiple stages with zero or more adaptive warm up stages followed by the main non-adaptive sampling stage. During the adaptive stage(s) parameters of the transition(s) are adaptively tuned based on the chain state and/or transition statistics.

The chains (including both adaptive and non-adaptive stages) may be run in parallel across multiple independent processes or sequentially. In all cases all chains use independent random draws.

Parameters:

n_warm_up_iter (int) – Number of adaptive warm up iterations per chain. Depending on the mici.stagers.Stager instance specified by the stager argument the warm up iterations may be split between one or more adaptive stages. If zero, only a single non-adaptive stage is used.
n_main_iter (int) – Number of iterations (samples to draw) per chain during main (non-adaptive) sampling stage.
init_states (Iterable[Union[ChainState, dict]]) – Initial chain states. Each entry can be either a ChainState object or a dictionary with entries specifying initial values for all state variables used by chain transition sample methods.
trace_funcs (Optional[Sequence[TraceFunction]]) – Sequence of functions which compute the variables to be recorded at each chain iteration (during only the main non-adaptive sampling stage if trace_warm_up == False), with each trace function passed the current state and returning a dictionary of scalar or array values corresponding to the variable(s) to be stored. The keys in the returned dictionaries are used to index the trace arrays in the returned traces dictionary. If a key appears in multiple dictionaries only the the value corresponding to the last trace function to return that key will be stored. If None or an empty sequence no variables are traced.
adapters (Optional[dict[str, Sequence[Adapter]]]) – Dictionary of sequences of mici.adapters.Adapter instances keyed by strings corresponding to the key of the transition in the transitions dictionary to apply the adapters to, to use to adaptatively set parameters of the transitions during the adaptive stages of the chains. Note that the adapter updates are applied in the order the adapters appear in the sequence and so if multiple adapters change the same parameter(s) the order will matter. If None or an empty sequence no adapters are used.
stager (Optional[Stager]) – Chain iteration stager object which controls the split of the chain iterations into the adaptive warm up and non-adaptive main stages. If set to None (the default) and all adapters specified by the adapters argument are of the fast type (i.e. their is_fast attribute is True) then a mici.stagers.WarmUpStager instance will be used corresponding to using a single adaptive warm up stage will all adapters active. If set to None and the adapters specified by the adapters argument are not all of the fast type, then a mici.stagers.WindowedWarmUpStager (with its default arguments) will be used, corresponding to using multiple adaptive warm up stages with only the fast-type adapters active in some - see documentation of mici.stagers.WarmUpStager for details.
n_process (Optional[int]) – Number of parallel processes to run chains over. If n_process=1 then chains will be run sequentially otherwise a multiprocessing.Pool object will be used to dynamically assign the chains across multiple processes. If set to None then the number of processes will be set to the output of os.cpu_count(). Default is n_process=1.
trace_warm_up (bool) – Whether to record chain traces and statistics during warm-up stage iterations (True) or only record traces and statistics in the iterations of the final non-adaptive stage (False, the default).
max_threads_per_process (Optional[int]) – If threadpoolctl is available this argument may be used to limit the maximum number of threads that can be used in thread pools used in libraries supported by threadpoolctl, which include BLAS and OpenMP implementations. This argument will only have an effect if n_process > 1 such that chains are being run on multiple processes and only if threadpoolctl is installed in the current Python environment. If set to None (the default) no limits are set.
force_memmap (bool) – Whether to force arrays used to store chain data to be memory-mapped to files on disk to avoid excessive system memory usage for long chains and/or large chain states. The chain data is written to .npy files in the directory specified by memmap_path (or a temporary directory if memmap_path is None). Chain data is always memory mapped when sampling chains in parallel on multiple processes.
memmap_path (Optional[str]) – Path to directory to write memory-mapped chain data to. If None (the default) and memory-mapping is enabled then a temporary directory will be created and the chain data written to files there, with the created files being deleted in this case once the last reference to them is closed.
monitor_stats (Optional[dict[str, list[str]]]) – String-keyed dictionary of lists of strings, with dictionary key the key of a Markov transition in the transitions dict passed to the the __init__() method and the corresponding list, the keys of statistics returned by the transition (as defined by the mici.transitions.Transition.statistics_type attribute of transition). The mean over samples computed so far of the statistics associated with any valid key-pairs will be monitored during sampling by printing as postfix to progress bar. If None no statistics are monitored.
display_progress (bool) – Whether to display a progress bar to track the completed chain sampling iterations. Default value is True, i.e. to display progress bar.
progress_bar_class (Optional[ProgressBar]) – Class or factory function for progress bar to use to show chain progress if enabled (display_progress=True). Defaults to mici.progressbars.SequenceProgressBar if None and display_progress=True.

Returns:

Named tuple (final_states, traces, statistics) corresponding to states of chains after final interatinos, dictionary of chain trace arrays and dictionary of chain statistics dictionaries.

Return type:

MCMCSampleChainsOutputs

property transitions: dict[str, Transition]#: Dictionary of transition kernels sampled from in each chain iteration.

class mici.samplers.RandomMetropolisHMC(system, integrator, rng, n_step_range, momentum_transition=None)[source]#

Bases: HamiltonianMonteCarlo

Random integration time HMC method with Metropolis sampling of new state.

In each transition a trajectory is generated by integrating the Hamiltonian dynamics from the current state in the current integration time direction for a random integer number of integrator steps sampled from the uniform distribution on an integer interval.

The state at the end of the trajectory with the integration direction negated (this ensuring the proposed move is an involution) is used as the proposal in a Metropolis acceptance step. The integration direction is then deterministically negated again irrespective of the accept decision, with the effect being that on acceptance the integration direction will be equal to its initial value and on rejection the integration direction will be the negation of its initial value.

The randomisation of the number of integration steps avoids the potential of the chain mixing poorly due to using an integration time close to the period of (near) periodic systems (Neal, 2011; Mackenzie, 1989).

References

Neal, R.M., 2011. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2(11), p.2.
Mackenzie, P.B., 1989. An improved hybrid Monte Carlo method. Physics Letters B, 226(3-4), pp.369-371.

Parameters:

system (System) – Hamiltonian system to be simulated.
rng (Generator) – Numpy random number generator.
integrator (Integrator) – Symplectic integrator to use to simulate dynamics in integration transition.
n_step_range (tuple[int, int]) – tuple (lower, upper) with two positive integer entries lower and upper (with upper > lower) specifying respectively the lower and upper bounds (inclusive) of integer interval to uniformly draw random number integrator steps to simulate in each integration transition.
momentum_transition (Optional[MomentumTransition]) – Markov transition kernel which leaves the conditional distribution on the momentum under the canonical distribution invariant, updating only the momentum component of the chain state. If set to None the momentum transition operator mici.transitions.IndependentMomentumTransition will be used, which independently samples the momentum from its conditional distribution.

property n_step_range: tuple[int, int]#: Interval to uniformly draw number of integrator steps from.

property rng: Generator#: NumPy random number generator object.

sample_chains(n_warm_up_iter, n_main_iter, init_states, **kwargs)#

Sample Markov chains from given initial states with adaptive warm up.

The default settings use a single (fast) DualAveragingStepSizeAdapter adapter instance which adapts the integrator step-size using a dual-averaging algorithm in a single adaptive stage.

The chains (including both adaptive and non-adaptive stages) may be run in parallel across multiple independent processes or sequentially. In all cases all chains use independent random draws.

Parameters:

n_warm_up_iter (int) – Number of adaptive warm up iterations per chain. Depending on the mici.stagers.Stager instance specified by the stager argument the warm up iterations may be split between one or more adaptive stages. If zero, only a single non-adaptive stage is used.
n_main_iter (int) – Number of iterations (samples to draw) per chain during main (non-adaptive) sampling stage.
init_states (Iterable[Union[ChainState, NDArray, dict]]) – Initial chain states. Each state can be either an array specifying the state position component or a mici.states.ChainState instance. If an array is passed or the mom attribute of the state is not set, a momentum component will be independently sampled from its conditional distribution. One chain will be run for each state in the iterable.

Keyword Arguments:

trace_funcs – Sequence of functions which compute the variables to be recorded at each chain iteration (during only the main non-adaptive sampling stage if trace_warm_up is False), with each trace function passed the current state and returning a dictionary of scalar or array values corresponding to the variable(s) to be stored. The keys in the returned dictionaries are used to index the trace arrays in the returned traces dictionary. If a key appears in multiple dictionaries only the the value corresponding to the last trace function to return that key will be stored. Default is to use a single function which records the position component of the state under the key "pos" and the Hamiltonian at the state under the key "hamiltonian".
adapters – Sequence of mici.adapters.Adapter instances to use to adaptatively set parameters of the integration transition such as the step size during the adaptive stages of the chains. Note that the adapter updates are applied in the order the adapters appear in the sequence and so if multiple adapters change the same parameter(s) the order will matter. If None or an empty sequence no adapters are used. Default is to use a single instance of mici.adapters.DualAveragingStepSizeAdapter with its default parameters.
stager – Chain iteration stager object which controls the split of the chain iterations into the adaptive warm up and non-adaptive main stages. If set to None (the default) and all adapters specified by the adapters argument are of the fast type (i.e. their is_fast attribute is True) then a mici.stagers.WarmUpStager instance will be used corresponding to using a single adaptive warm up stage will all adapters active. If set to None and the adapters specified by the adapters argument are not all of the fast type, then a mici.stagers.WindowedWarmUpStager (with its default arguments) will be used, corresponding to using multiple adaptive warm up stages with only the fast-type adapters active in some - see documentation of mici.stagers.WarmUpStager for details.
n_process – Number of parallel processes to run chains over. If n_process=1 then chains will be run sequentially otherwise a multiprocessing.Pool object will be used to dynamically assign the chains across multiple processes. If set to None then the number of processes will be set to the output of os.cpu_count(). Default is n_process=1.
trace_warm_up – Whether to record chain traces and statistics during warm-up stage iterations (True) or only record traces and statistics in the iterations of the final non-adaptive stage (False, the default).
max_threads_per_process – If threadpoolctl is available this argument may be used to limit the maximum number of threads that can be used in thread pools used in libraries supported by threadpoolctl, which include BLAS and OpenMP implementations. This argument will only have an effect if n_process > 1 such that chains are being run on multiple processes and only if threadpoolctl is installed in the current Python environment. If set to None (the default) no limits are set.
force_memmap – Whether to force arrays used to store chain data to be memory-mapped to files on disk to avoid excessive system memory usage for long chains and/or large chain states. The chain data is written to .npy files in the directory specified by memmap_path (or a temporary directory if memmap_path is None). Chain data is always memory mapped when sampling chains in parallel on multiple processes.
memmap_path – Path to directory to write memory-mapped chain data to. If None (the default) and memory-mapping is enabled then a temporary directory will be created and the chain data written to files there, with the created files being deleted in this case once the last reference to them is closed.
monitor_stats – Sequence of string keys of (integration) transition statistics to monitor mean of over samples computed so far during sampling by printing as postfix to progress bar. Default is to print only the mean accept_stat (acceptance statistic).
display_progress – Whether to display a progress bar to track the completed chain sampling iterations. Default value is True, i.e. to display progress bar.
progress_bar_class – Class or factory function for progress bar to use to show chain progress if enabled (display_progress=True). Defaults to mici.progressbars.SequenceProgressBar.

Returns:

Named tuple (final_states, traces, statistics) corresponding to states of chains after final iteration, dictionary of chain trace arrays and dictionary of chain statistics dictionaries.

Return type:

HMCSampleChainsOutputs

property system: System#: Hamiltonian system corresponding to joint distribution on augmented space.

property transitions: dict[str, Transition]#: Dictionary of transition kernels sampled from in each chain iteration.

class mici.samplers.StaticMetropolisHMC(system, integrator, rng, n_step, momentum_transition=None)[source]#

Bases: HamiltonianMonteCarlo

Static integration time HMC method with Metropolis sampling.

In each transition a trajectory is generated by integrating the Hamiltonian dynamics from the current state in the current integration time direction for a fixed integer number of integrator steps.

This is original proposed Hybrid Monte Carlo (often now instead termed Hamiltonian Monte Carlo) algorithm (Duane et al., 1987; Neal, 2011).

References

Duane, S., Kennedy, A.D., Pendleton, B.J. and Roweth, D., 1987. Hybrid Monte Carlo. Physics letters B, 195(2), pp.216-222.
Neal, R.M., 2011. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2(11), p.2.

Parameters:

system (System) – Hamiltonian system to be simulated.
rng (Generator) – Numpy random number generator.
integrator (Integrator) – Symplectic integrator to use to simulate dynamics in integration transition.
n_step (int) – Number of integrator steps to simulate in each integration transition.
momentum_transition (Optional[MomentumTransition]) – Markov transition kernel which leaves the conditional distribution on the momentum under the canonical distribution invariant, updating only the momentum component of the chain state. If set to None the momentum transition operator mici.transitions.IndependentMomentumTransition will be used, which independently samples the momentum from its conditional distribution.

property n_step: int#: Number of integrator steps per integrator transition.

property rng: Generator#: NumPy random number generator object.

sample_chains(n_warm_up_iter, n_main_iter, init_states, **kwargs)#

Sample Markov chains from given initial states with adaptive warm up.

The default settings use a single (fast) DualAveragingStepSizeAdapter adapter instance which adapts the integrator step-size using a dual-averaging algorithm in a single adaptive stage.

The chains (including both adaptive and non-adaptive stages) may be run in parallel across multiple independent processes or sequentially. In all cases all chains use independent random draws.

Parameters:

n_warm_up_iter (int) – Number of adaptive warm up iterations per chain. Depending on the mici.stagers.Stager instance specified by the stager argument the warm up iterations may be split between one or more adaptive stages. If zero, only a single non-adaptive stage is used.
n_main_iter (int) – Number of iterations (samples to draw) per chain during main (non-adaptive) sampling stage.
init_states (Iterable[Union[ChainState, NDArray, dict]]) – Initial chain states. Each state can be either an array specifying the state position component or a mici.states.ChainState instance. If an array is passed or the mom attribute of the state is not set, a momentum component will be independently sampled from its conditional distribution. One chain will be run for each state in the iterable.

Keyword Arguments:

trace_funcs – Sequence of functions which compute the variables to be recorded at each chain iteration (during only the main non-adaptive sampling stage if trace_warm_up is False), with each trace function passed the current state and returning a dictionary of scalar or array values corresponding to the variable(s) to be stored. The keys in the returned dictionaries are used to index the trace arrays in the returned traces dictionary. If a key appears in multiple dictionaries only the the value corresponding to the last trace function to return that key will be stored. Default is to use a single function which records the position component of the state under the key "pos" and the Hamiltonian at the state under the key "hamiltonian".
adapters – Sequence of mici.adapters.Adapter instances to use to adaptatively set parameters of the integration transition such as the step size during the adaptive stages of the chains. Note that the adapter updates are applied in the order the adapters appear in the sequence and so if multiple adapters change the same parameter(s) the order will matter. If None or an empty sequence no adapters are used. Default is to use a single instance of mici.adapters.DualAveragingStepSizeAdapter with its default parameters.
stager – Chain iteration stager object which controls the split of the chain iterations into the adaptive warm up and non-adaptive main stages. If set to None (the default) and all adapters specified by the adapters argument are of the fast type (i.e. their is_fast attribute is True) then a mici.stagers.WarmUpStager instance will be used corresponding to using a single adaptive warm up stage will all adapters active. If set to None and the adapters specified by the adapters argument are not all of the fast type, then a mici.stagers.WindowedWarmUpStager (with its default arguments) will be used, corresponding to using multiple adaptive warm up stages with only the fast-type adapters active in some - see documentation of mici.stagers.WarmUpStager for details.
n_process – Number of parallel processes to run chains over. If n_process=1 then chains will be run sequentially otherwise a multiprocessing.Pool object will be used to dynamically assign the chains across multiple processes. If set to None then the number of processes will be set to the output of os.cpu_count(). Default is n_process=1.
trace_warm_up – Whether to record chain traces and statistics during warm-up stage iterations (True) or only record traces and statistics in the iterations of the final non-adaptive stage (False, the default).
max_threads_per_process – If threadpoolctl is available this argument may be used to limit the maximum number of threads that can be used in thread pools used in libraries supported by threadpoolctl, which include BLAS and OpenMP implementations. This argument will only have an effect if n_process > 1 such that chains are being run on multiple processes and only if threadpoolctl is installed in the current Python environment. If set to None (the default) no limits are set.
force_memmap – Whether to force arrays used to store chain data to be memory-mapped to files on disk to avoid excessive system memory usage for long chains and/or large chain states. The chain data is written to .npy files in the directory specified by memmap_path (or a temporary directory if memmap_path is None). Chain data is always memory mapped when sampling chains in parallel on multiple processes.
memmap_path – Path to directory to write memory-mapped chain data to. If None (the default) and memory-mapping is enabled then a temporary directory will be created and the chain data written to files there, with the created files being deleted in this case once the last reference to them is closed.
monitor_stats – Sequence of string keys of (integration) transition statistics to monitor mean of over samples computed so far during sampling by printing as postfix to progress bar. Default is to print only the mean accept_stat (acceptance statistic).
display_progress – Whether to display a progress bar to track the completed chain sampling iterations. Default value is True, i.e. to display progress bar.
progress_bar_class – Class or factory function for progress bar to use to show chain progress if enabled (display_progress=True). Defaults to mici.progressbars.SequenceProgressBar.

Returns:

Named tuple (final_states, traces, statistics) corresponding to states of chains after final iteration, dictionary of chain trace arrays and dictionary of chain statistics dictionaries.

Return type:

HMCSampleChainsOutputs

property system: System#: Hamiltonian system corresponding to joint distribution on augmented space.

property transitions: dict[str, Transition]#: Dictionary of transition kernels sampled from in each chain iteration.