mot.sample package

Submodules

mot.sample.amwg module

class mot.sample.amwg.AdaptiveMetropolisWithinGibbs(ll_func, log_prior_func, x0, proposal_stds, target_acceptance_rate=0.44, batch_size=50, damping_factor=1, min_val=1e-15, max_val=1000.0, **kwargs)[source]

Bases: mot.sample.base.AbstractRWMSampler

An implementation of the Adaptive Metropolis-Within-Gibbs (AMWG) MCMC algorithm [1].

This scales the proposal parameter (typically the std) such that it oscillates towards the chosen acceptance rate. We implement the delta function (see [1]) as: \(\delta(n) = \sqrt{1 / (d*n)}\). Where n is the current batch index and d is the damping factor. As an example, with a damping factor of 500, delta reaches a scaling of 0.01 in 20 batches. At a batch size of 50 that would amount to 1000 samples.

In principal, AMWG can be used as a final MCMC algorithm if and only if the adaption is effectively zero. That is when delta gets close enough to zero to no longer influence the proposals.

Parameters:
  • ll_func (mot.lib.cl_function.CLFunction) – The log-likelihood function. See parent docs.
  • log_prior_func (mot.lib.cl_function.CLFunction) – The log-prior function. See parent docs.
  • x0 (ndarray) – the starting positions for the sampler. Should be a two dimensional matrix with for every modeling instance (first dimension) and every parameter (second dimension) a value.
  • proposal_stds (ndarray) – for every parameter and every modeling instance an initial proposal std.
  • target_acceptance_rate (float) – the target acceptance rate between 0 and 1.
  • batch_size (int) – the size of the batches in between which we update the parameters
  • damping_factor (int) – how fast the adaptation moves to zero
  • min_val (float) – the minimum value the standard deviation can take
  • max_val (float) – the maximum value the standard deviation can take

References

[1] Roberts GO, Rosenthal JS. Examples of adaptive MCMC. J Comput Graph Stat. 2009;18(2):349-367.
doi:10.1198/jcgs.2009.06134.

mot.sample.base module

class mot.sample.base.AbstractRWMSampler(ll_func, log_prior_func, x0, proposal_stds, use_random_scan=False, finalize_proposal_func=None, **kwargs)[source]

Bases: mot.sample.base.AbstractSampler

An abstract basis for Random Walk Metropolis (RWM) samplers.

Random Walk Metropolis (RWM) samplers require for every parameter and every modeling instance an proposal standard deviation, used in the random walk.

Parameters:
  • ll_func (mot.lib.cl_function.CLFunction) – The log-likelihood function.
  • log_prior_func (mot.lib.cl_function.CLFunction) – The log-prior function.
  • x0 (ndarray) – the starting positions for the sampler. Should be a two dimensional matrix with for every modeling instance (first dimension) and every parameter (second dimension) a value.
  • proposal_stds (ndarray) – for every parameter and every modeling instance an initial proposal std.
  • use_random_scan (boolean) – if we iterate over the parameters in a random order or in a linear order at every sample iteration. By default we apply a system scan starting from the first dimension to the last. With a random scan we randomize the indices every iteration.
  • finalize_proposal_func (mot.lib.cl_function.CLFunction) –

    a CL function to finalize every proposal by the sampling routine. This allows the model to change a proposal before computing the prior or likelihood probabilities. If None, we will not use this callback.

    As an example, suppose you are sampling a polar coordinate :math:` heta` defined on \([0, 2\pi]\) with a random walk Metropolis proposal distribution. This distribution might propose positions outside of the range of :math:` heta`. Of course the model function could deal with that by taking the modulus of the input, but then you have to post-process the chain with the same transformation. Instead, this function allows changing the proposal before it is put into the model and before it will be stored in the chain.

    This function should return a proposal that is equivalent (but not necessarily equal) to the provided proposal.

    Signature:

    void <func_name>(void* data, local mot_float_type* x);
class mot.sample.base.AbstractSampler(ll_func, log_prior_func, x0, data=None, cl_runtime_info=None, **kwargs)[source]

Bases: object

Abstract base class for sample routines.

Sampling routines implementing this interface should be stateful objects that, for the given likelihood and prior, keep track of the sample state over multiple calls to sample().

Parameters:
  • ll_func (mot.lib.cl_function.CLFunction) –

    The log-likelihood function. A CL function with the signature:

    double <func_name>(local const mot_float_type* const x, void* data);
  • log_prior_func (mot.lib.cl_function.CLFunction) –

    The log-prior function. A CL function with the signature:

    mot_float_type <func_name>(local const mot_float_type* const x, void* data);
  • x0 (ndarray) – the starting positions for the sampler. Should be a two dimensional matrix with for every modeling instance (first dimension) and every parameter (second dimension) a value.
  • data (mot.lib.kernel_data.KernelData) – the user provided data for the void* data pointer.
sample(nmr_samples, burnin=0, thinning=1)[source]

Take additional samples from the given likelihood and prior, using this sampler.

This method can be called multiple times in which the sample state is stored in between.

Parameters:
  • nmr_samples (int) – the number of samples to return
  • burnin (int) – the number of samples to discard before returning samples
  • thinning (int) – how many sample we wait before storing a new one. This will draw extra samples such that the total number of samples generated is nmr_samples * (thinning) and the number of samples stored is nmr_samples. If set to one or lower we store every sample after the burn in.
Returns:

the sample output object
Return type:

SamplingOutput
set_cl_runtime_info(cl_runtime_info)[source]

Update the CL runtime information.

Parameters:cl_runtime_info (mot.configuration.CLRuntimeInfo) – the new runtime information
class mot.sample.base.SamplingOutput[source]

Bases: object

get_log_likelihoods()[source]

Get per set of sampled parameters the log likelihood value associated with that set of parameters.

Returns:
the log likelihood values, a (d, n) array with for d problems and n samples the log likelihood
value.
Return type:ndarray
get_log_priors()[source]

Get per set of sampled parameters the log prior value associated with that set of parameters.

Returns:the log prior values, a (d, n) array with for d problems and n samples the prior value.
Return type:ndarray
get_samples()[source]

Get the matrix containing the sample results.

Returns:the sampled parameter maps, a (d, p, n) array with for d problems and p parameters n samples.
Return type:ndarray
class mot.sample.base.SimpleSampleOutput(samples, log_likelihoods, log_priors)[source]

Bases: mot.sample.base.SamplingOutput

Simple storage container for the sample output

get_log_likelihoods()[source]

Get per set of sampled parameters the log likelihood value associated with that set of parameters.

Returns:
the log likelihood values, a (d, n) array with for d problems and n samples the log likelihood
value.
Return type:ndarray
get_log_priors()[source]

Get per set of sampled parameters the log prior value associated with that set of parameters.

Returns:the log prior values, a (d, n) array with for d problems and n samples the prior value.
Return type:ndarray
get_samples()[source]

Get the matrix containing the sample results.

Returns:the sampled parameter maps, a (d, p, n) array with for d problems and p parameters n samples.
Return type:ndarray

mot.sample.mwg module

class mot.sample.mwg.MetropolisWithinGibbs(ll_func, log_prior_func, x0, proposal_stds, **kwargs)[source]

Bases: mot.sample.base.AbstractRWMSampler

An implementation of the Metropolis-Within-Gibbs MCMC algorithm [1].

This does not scale the proposal standard deviations during sampling.

Parameters:
  • ll_func (mot.lib.cl_function.CLFunction) – The log-likelihood function. See parent docs.
  • log_prior_func (mot.lib.cl_function.CLFunction) – The log-prior function. See parent docs.
  • x0 (ndarray) – the starting positions for the sampler. Should be a two dimensional matrix with for every modeling instance (first dimension) and every parameter (second dimension) a value.
  • proposal_stds (ndarray) – for every parameter and every modeling instance an initial proposal std.

References

[1] 1. van Ravenzwaaij D, Cassey P, Brown SD. A simple introduction to Markov Chain Monte–Carlo sampling.
Psychon Bull Rev. 2016:1-12. doi:10.3758/s13423-016-1015-8.

mot.sample.scam module

class mot.sample.scam.SingleComponentAdaptiveMetropolis(ll_func, log_prior_func, x0, proposal_stds, waiting_period=100, scaling_factor=2.4, epsilon=1e-05, **kwargs)[source]

Bases: mot.sample.base.AbstractRWMSampler

An implementation of the Single Component Adaptive Metropolis (SCAM) MCMC algorithm [1].

The SCAM method works by adapting the proposal standard deviation to the empirical standard deviation of the component’s marginal distribution. That is, the standard deviation \(\sigma_i^{(t)}\) for the proposal distribution of the \(i\) th component at time \(t\) is given by:

\[\begin{split}\sigma_i^{(t)} = \begin{cases} \sigma_i^{(0)}, & t \leq t_s \\ 2.4 * \sqrt{\mathrm{Var}(\mathbf{X}^{(0)}_i, \ldots, \mathbf{X}^{(t-1)}_i)} + 1\cdot \epsilon, & t > t_s \end{cases}\end{split}\]

where \(t_s\) denotes the iteration after which the adaptation starts (we use \(t_s = 100\)). A small constant is necessary to prevent the standard deviation from shrinking to zero. This adaptation algorithm has been proven to retain ergodicity, meaning it is guaranteed to converge to the right stationary distribution [1].

Parameters:
  • ll_func (mot.lib.cl_function.CLFunction) – The log-likelihood function. See parent docs.
  • log_prior_func (mot.lib.cl_function.CLFunction) – The log-prior function. See parent docs.
  • x0 (ndarray) – the starting positions for the sampler. Should be a two dimensional matrix with for every modeling instance (first dimension) and every parameter (second dimension) a value.
  • proposal_stds (ndarray) – for every parameter and every modeling instance an initial proposal std.
  • waiting_period (int) – only start updating the proposal std. after this many draws.
  • scaling_factor (float) – the scaling factor to use (the parameter s in the paper referenced).
  • epsilon (float or ndarray) – small number to prevent the std. from collapsing to zero. Can either be one value for all parameters, or one value per parameter.

References

[1] Haario, H., Saksman, E., & Tamminen, J. (2005). Componentwise adaptation for high dimensional MCMC.
Computational Statistics, 20(2), 265-273. https://doi.org/10.1007/BF02789703

Module contents