samplers

Users typically do not need to access samplers directly. Instead, they inform calibrators which sampler should be used and provide the calibrators with the full set of configuration values required by their sampler. This interface information is provided so that users can determine which sampler they would like to use and how to configure it for their needs.

surmise.utilitiesmethods.sample_with_metropolis_hastings(logpost_func, draw_func, scipy_stats_rng, specification)

Metropolis Hastings Markov chain Monte Carlo sampling method.

Parameters:

  • logpost_func (function) – function that returns the log of the posterior for a given theta provided as a \(1 \times p\) 2D NumPy row vector.

  • draw_func (function) – function that accepts the number of desired random draws needed for initializing the sample process and returns a 2D NumPy array of draws with each row being a different draw.

  • scipy_stats_rng – scipy.stats-compatible pseudorandom number generator that the sampler should use for all random draws performed by the sampler. The sampling process produces identical results if it is repeated with the same RNG setup.

  • specification (dict) –

    The full set of sampler configuration values

    • ”theta0” - None or the initial theta to use to start the sampling process. If None, then 100 candidates are drawn using draw_func and the first with finite log posterior is used. An error is raised if no candidates has a finite log posterior.

    • ”nSamples” - total number of samples to acquire after the burn-in period.

    • ”nBurnSamples” - total number of samples to acquire during the burn-in period.

    • ”stepType” - a multivariate uniform step proposal distribution centered on zero is used if “uniform” is provided; a zero-mean multivariate normal step proposal distribution, if “normal” is provided.

    • ”stepParam” - None or the lengthscales that characterize the step proposal distribution.

      • widths of uniform distribution if stepType is “uniform”

      • standard deviations of multivariate normal distribution if stepType is “normal”

      Note that for “normal” the covariances are all set to zero. If None, then the lengthscale is set to the standard deviations of nBurnSamples random draws from draw_func.

    • ”verbose” - log setup and sampling progress information if True.

Returns:

sampler_info – Summary of the sampling process

  • ”theta” - 2D NumPy array whose rows are the accepted theta samples provided in the order in which they were determined. This does not include theta determined during the burn-in period.

  • ”lpostlist” - 1D NumPy array of log posterior values obtained at all candidate theta, including those rejected by the sampling process, provided in the order of evalution. This includes the values obtained during the burn-in period.

  • ”acc_rate” - final acceptance rate of the sampling process derived from the determination of only the final nSamples theta

Return type:

dict

Todo

  • Once the sampler arguments have been separated out from all other arguments in higher-level code, the samplers should confirm that they are passed values for all arguments and no more.

surmise.utilitiesmethods.sample_with_LMC(logpost_func, draw_func, scipy_stats_rng, specification)

Note

This sampler is currently considered as research-grade and is not officially offered by surmise.

Metropolis-adjusted Langevin algorithm or Langevin Monte Carlo (LMC), which seeks to propose the next iterates by leveraging gradient information at the current iterate. The proposal has the form

\[\theta^{k+1} = \theta^k - \nabla g(\theta^k) \Delta t + \sqrt{2\Delta t} Z,\]

where \(\Delta t\) is a time stepsize, and \(Z\) is an independently and identically drawn sample from the standard Gaussian normal of the appropriate dimension. The proposal is then accepted or rejected by the typical Metropolis-Hastings step, i.e. accept with probability

\[\alpha = \min\left\{1, \frac{\pi(\tilde{\theta}^{k+1})q(\theta^k \mid \tilde{\theta}^{k+1})}{\pi(\theta^{k})q(\tilde{\theta}^{k+1} \mid \theta^k)}\right\},\]

where \(\pi(\cdot)\) is the posterior distribution, \(q(\cdot \mid \cdot)\) is the proposal distribution, and \(\theta^k, \tilde{\theta}^{k+1}\) are the current and the proposed point respectively.

Langevin Monte Carlo has shown strengths in increasing the acceptance rate, compared to the typical Metropolis-Hastings algorithm (Roberts and Rosenthal, 1998). However, its significant drawback lies in its poor scaling due to the computation for the gradient at the current iterate.

Refer to G. O. Roberts and J. S. Rosenthal. Optimal scaling of discrete approximations to langevin diffusions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255-268, 1998.

Parameters:

  • logpostfunc (function) – A function that returns the log of the posterior densities at each of \(m\) theta points provided in an \(m \times p\) NumPy array. If gradients are not computed, logpostfunc should return a length \(m\) NumPy array of log posterior values computed at the given theta points. If gradients are computed, logpostfunc should return a tuple whose first element is the array of log posterior values as described above; the second, an \(m \times p\) NumPy array of gradients of the log posterior at those same theta points.

  • draw_func (function) – function that accepts the number of desired random draws needed for initializing the sample process and returns a 2D NumPy array of draws with each row being a different draw.

  • scipy_stats_rng – scipy.stats-compatible pseudorandom number generator that the sampler should use for all random draws performed by the sampler. The sampling process produces identical results if it is repeated with the same RNG setup.

  • specification (dict) –

    The full set of sampler configuration values

    • ”theta0” - None or an \(m \times p\) array of initial thetas to use to start the sampling process. If None, then the process is intialized with 1000 random draws from draw_func.

    • ”nSamples” - total number of samples from the posterior.

    • ”verbose” - log setup and sampling progress information if True.

Returns:

sampler_info – Summary of the sampling process

  • ”theta” - an nSamples \(\times p\) array of unsorted samples from the posterior

  • ”logpost” - an nSamples-element array of log posterior values associated with the elements of theta

Return type:

dict

surmise.utilitiesmethods.sample_with_PTLMC(logpost_func, draw_func, scipy_stats_rng, specification)

Parallel-Tempering Ensemble Markov chain Monte Carlo sampling method based on Langevin Monte Carlo (See surmise.utilitiesmethods.sample_with_LMC()).

Parameters:

  • logpost_func (function) – A function that returns the log of the posterior densities at each of \(m\) theta points provided in an \(m \times p\) NumPy array. If gradients are not computed, logpost_func should return a length \(m\) NumPy array of log posterior values computed at the given theta points. If gradients are computed, logpost_func should return a tuple whose first element is the array of log posterior values as described above; the second, an \(m \times p\) NumPy array of gradients of the log posterior at those same theta points.

  • draw_func (function) – function that accepts the number of desired random draws needed for initializing the sample process and returns a 2D NumPy array of draws with each row being a different draw.

  • scipy_stats_rng – scipy.stats-compatible pseudorandom number generator that the sampler should use for all random draws performed by the sampler. The sampling process produces identical results if it is repeated with the same RNG setup.

  • specification (dict) –

    The full set of sampler configuration values

    • ”theta0” - None or an \(m \times p\) array of initial thetas to use to start the sampling process. If None or too few initial theta are provided, the process is intialized with 1000 random draws from draw_func.

    • ”nSamples” - positive integer that controls how many samples are returned in theta

    • ”nChains” - positive integer that controls how many chains of fixed temperature to run simultaneously.

    • ”samplesPerChain” - positive integer that controls how many samples should be made for each chain.

    • ”nTemperatures” - positive integer that controls how many chains of varying temperature to run simultaneously.

    • ”maxTemperature” - number greater than 1 that gives the maximum temperature used in parallel tempering.

    • ”verbose” - log setup and sampling progress information if True.

Returns:

sampler_info – Summary of the sampling process

Todo

revisit when flattening is addressed.

  • ”theta” - an nSamples \(\times p\) array: the first nSamples entries from the flattened samples [(chain 1 … chain 2 … chain nChains)]

  • ”theta_from_chain” - an nChains \(\times\) samplesPerChains \(\times p\) array of unflattened, unsorted samples.

  • ”logpost” - an nSamples-element array of log posterior values associated with the elements of theta

Return type:

dict