pymcmcstat.samplers package

pymcmcstat.samplers.Adaptation module

Created on Thu Jan 18 11:14:11 2018

@author: prmiles

class pymcmcstat.samplers.Adaptation.Adaptation[source]

Adaptive Metropolis (AM) algorithm based on [haario2001adaptive].

Attributes:
run_adaptation(covariance, options, isimu, iiadapt, rejected, chain, chainind, u, npar, alpha)[source]

Run adaptation step

Args:
  • covariance (CovarianceProcedures): Covariance methods and variables
  • options (SimulationOptions): Options for MCMC simulation
  • isimu (int): Simulation counter
  • iiadapt (int): Adaptation counter
  • rejected (dict): Rejection counter
  • chain (ndarray): Sampling chain
  • chainind (ind): Relative point in chain
  • u (ndarray): Latest random sample points
  • npar (int): Number of parameters being sampled
  • alpha (float): Latest Likelihood evaluation
Returns:
pymcmcstat.samplers.Adaptation.adjust_cov_matrix(upcov, R, npar, qcov_adjust, qcov_scale, rejected, iiadapt, verbosity)[source]

Adjust covariance matrix if found to be singular.

Args:
  • upcov (ndarray): Parameter covariance matrix.
  • R (ndarray): Cholesky decomposition of covariance matrix.
  • npar (int): Number of parameters.
  • qcov_adjust (float): Covariance adjustment factor.
  • qcov_scale (float): Scale parameter
  • rejected (dict): Rejection counters.
  • iiadapt (int): Adaptation counter.
  • verbosity (int): Verbosity of display output.
Returns:
  • R (ndarray): Cholesky decomposition of covariance matrix.
pymcmcstat.samplers.Adaptation.below_burnin_threshold(rejected, iiadapt, R, burnin_scale, verbosity)[source]

Update Cholesky Matrix using below burnin thershold

Args:
  • rejected (dict): Rejection counters.
  • iiadapt (int): Adaptation counter
  • R (ndarray): Cholesky decomposition of covariance matrix.
  • burnin_scale (float): Scale for burnin.
  • verbosity (int): Verbosity of display output.
Returns:
  • R (ndarray): Scaled Cholesky matrix.
pymcmcstat.samplers.Adaptation.check_for_singular_cov_matrix(upcov, R, npar, qcov_adjust, qcov_scale, rejected, iiadapt, verbosity)[source]

Check if singular covariance matrix

Args:
  • upcov (ndarray): Parameter covariance matrix.
  • R (ndarray): Cholesky decomposition of covariance matrix.
  • npar (int): Number of parameters.
  • qcov_adjust (float): Covariance adjustment factor.
  • qcov_scale (float): Scale parameter
  • rejected (dict): Rejection counters.
  • iiadapt (int): Adaptation counter.
  • verbosity (int): Verbosity of display output.
Returns:
  • R (ndarray): Adjusted Cholesky decomposition of covariance matrix.
pymcmcstat.samplers.Adaptation.cholupdate(R, x)[source]

Update Cholesky decomposition

Args:
  • R (ndarray): Weighted Cholesky decomposition
  • x (ndarray): Weighted sum based on local chain update
Returns:
  • R1 (ndarray): Updated Cholesky decomposition
pymcmcstat.samplers.Adaptation.initialize_covariance_mean_sum(x, w)[source]

Initialize covariance chain, local mean, local sum

Args:
Returns:
  • xcov (ndarray): Initial covariance matrix
  • xmean (ndarray): Initial mean chain values
  • wsum (ndarray): Initial weighted sum
pymcmcstat.samplers.Adaptation.is_semi_pos_def_chol(x)[source]

Check if matrix is semi-positive definite using Cholesky Decomposition

Args:
Returns:
  • Boolean
  • c (ndarray): Cholesky decomposition (upper triangular form) or None
pymcmcstat.samplers.Adaptation.scale_cholesky_decomposition(Ra, qcov_scale)[source]

Scale Cholesky decomposition

Args:
  • R (ndarray): Cholesky decomposition of covariance matrix.
  • qcov_scale (float): Scale factor
Returns:
  • R (ndarray): Scaled Cholesky decomposition of covariance matrix.
pymcmcstat.samplers.Adaptation.setup_cholupdate(R, oldwsum, wsum, oldmean, xi)[source]

Setup input arguments to the Cholesky update function

Args:
  • R (ndarray): Previous Cholesky decomposition matrix
  • oldwsum (ndarray): Previous weighted sum
  • w (ndarray): Current Weights
  • oldmean (ndarray): Previous mean chain values
  • xi (ndarray): Row of chain segment
Returns:
  • Rin (ndarray): Scaled Cholesky decomposition matrix
  • xin (ndarray): Updated mean chain values for Cholesky function
pymcmcstat.samplers.Adaptation.setup_w_R(w, oldR, n)[source]

Setup weights and Cholesky matrix

Args:
  • x (ndarray): Chain segment
  • w (ndarray): Weights
  • oldcov (ndarray or None): Previous covariance matrix
  • oldmean (ndarray): Previous mean chain values
  • oldwsum (ndarray): Previous weighted sum
  • oldR (ndarray): Previous Cholesky decomposition matrix
Returns:
  • w (ndarray): Weights
  • R (ndarray): Previous Cholesky decomposition matrix
pymcmcstat.samplers.Adaptation.unpack_covariance_settings(covariance)[source]

Unpack covariance settings

Args:
Returns:
  • last_index_since_adaptation (int): Last index since adaptation occured.
  • R (ndarray): Cholesky decomposition of covariance matrix.
  • oldcovchain (ndarray): Covariance matrix history.
  • oldmeanchain (ndarray): Current mean chain values.
  • oldwsum (ndarray): Weights
  • no_adapt_index (numpy.ndarray): Indices of parameters not being adapted.
  • qcov_scale (float): Scale parameter
  • qcov (ndarray): Covariance matrix
pymcmcstat.samplers.Adaptation.unpack_simulation_options(options)[source]

Unpack simulation options

Args:
Returns:
  • burnintime (int):
  • burnin_scale (float): Scale for burnin.
  • ntry (int): Number of tries to take before rejection. Default is method dependent.
  • drscale (ndarray): Reduced scale for sampling in DR algorithm. Default is [5,4,3].
  • alphatarget (float): Acceptance ratio target.
  • etaparam (float):
  • qcov_adjust (float): Adjustment scale for covariance matrix.
  • doram (int): Flag to perform 'ram' algorithm (Obsolete).
  • verbosity (int): Verbosity of display output.
pymcmcstat.samplers.Adaptation.update_cov_from_covchain(covchain, qcov, no_adapt_index)[source]

Update covariance matrix from covariance matrix chain

Args:
  • covchain (ndarray): Covariance matrix history.
  • qcov (ndarray): Parameter covariance matrix
  • no_adapt_index (numpy.ndarray): Indices of parameters not being adapted.
Returns:
  • upcov (ndarray): Updated covariance matrix
pymcmcstat.samplers.Adaptation.update_cov_via_ram(u, isimu, etaparam, npar, alphatarget, alpha, R)[source]

Update covariance matrix via RAM

Args:
  • u (ndarray): Latest random sample points
  • isimu (int): Simulation counter
  • alphatarget (float): Acceptance ratio target.
  • npar (int): Number of parameters.
  • etaparam (float):
  • alpha (float): Latest Likelihood evaluation
  • R (ndarray): Cholesky decomposition of covariance matrix.
Returns:
  • upcov (ndarray): Updated parameter covariance matrix.
pymcmcstat.samplers.Adaptation.update_covariance_mean_sum(x, w, oldcov, oldmean, oldwsum, oldR=None)[source]

Update covariance chain, local mean, local sum

Args:
  • x (ndarray): Chain segment
  • w (ndarray): Weights
  • oldcov (ndarray or None): Previous covariance matrix
  • oldmean (ndarray): Previous mean chain values
  • oldwsum (ndarray): Previous weighted sum
  • oldR (ndarray): Previous Cholesky decomposition matrix
Returns:
  • xcov (ndarray): Updated covariance matrix
  • xmean (ndarray): Updated mean chain values
  • wsum (ndarray): Updated weighted sum
pymcmcstat.samplers.Adaptation.update_delayed_rejection(R, npar, ntry, drscale)[source]

Update Cholesky/Inverse matrix for Delayed Rejection

Args:
  • R (ndarray): Cholesky decomposition of covariance matrix.
  • npar (int): Number of parameters.
  • ntry (int): Number of tries to take before rejection. Default is method dependent.
  • drscale (ndarray): Reduced scale for sampling in DR algorithm. Default is [5,4,3].
Returns:
  • RDR (list): List of Cholesky decomposition of covariance matrices for each stage of DR.
  • InvR (list): List of Inverse Cholesky decomposition of covariance matrices for each stage of DR.

pymcmcstat.samplers.DelayedRejection module

Created on Thu Jan 18 10:42:07 2018

@author: prmiles

class pymcmcstat.samplers.DelayedRejection.DelayedRejection[source]

Delayed Rejection (DR) algorithm based on [haario2006dram].

Attributes:
alphafun(trypath, invR)[source]

Calculate likelihood according to DR

Args:
  • trypath (list): Sequence of DR steps
  • invR (ndarray): Inverse Cholesky decomposition matrix
Returns:
  • alpha (float): Result of likelihood function according to delayed rejection
classmethod initialize_next_metropolis_step(npar, old_theta, sigma2, RDR)[source]

Take metropolis step according to DR

Args:
  • npar (int): Number of parameters
  • old_theta (ndarray): q^{k-1}
  • sigma2 (float): Observation error variance
  • RDR (ndarray): Cholesky decomposition of parameter covariance matrix for DR steps
  • itry (int): DR step counter
Returns:
  • next_set (ParameterSet): New proposal set
  • u (numpy.ndarray): Numbers sampled from standard normal distributions (u.shape = (1,npar))
run_delayed_rejection(old_set, new_set, RDR, ntry, parameters, invR, sosobj, priorobj)[source]

Perform delayed rejection step - occurs in standard metropolis is not accepted.

Args:
Returns:
  • accept (int): 0 - reject, 1 - accept
  • out_set (ParameterSet): If accept == 1, then latest DR set; Else, q^k=q^{k-1}
  • outbound (int): 1 - rejected due to sampling outside of parameter bounds
pymcmcstat.samplers.DelayedRejection.extract_state_elements(iq, stage, trypath)[source]

Extract elements from tried paths.

Args:
  • iq (int): Stage number.
  • stage (int): Number of stages - 2
  • trypath (list): Sequence of DR steps
pymcmcstat.samplers.DelayedRejection.log_posterior_ratio(x1, x2)[source]

Calculate the logarithm of the posterior ratio.

Args:
Returns:
  • zq (float): Logarithm of posterior ratio.
pymcmcstat.samplers.DelayedRejection.nth_stage_log_proposal_ratio(iq, trypath, invR)[source]

Gaussian nth stage log proposal ratio.

Logarithm of q_i(y_n,...,y_{n-j})/q_i(x,y_1,...,y_j)

Args:
  • iq (int): Stage number.
  • trypath (list): Sequence of DR steps
  • invR (ndarray): Inverse Cholesky decomposition matrix
Returns:
  • zq (float): Logarithm of Gaussian nth stage proposal ratio.
pymcmcstat.samplers.DelayedRejection.update_set_based_on_acceptance(accept, old_set, next_set)[source]

Define output set based on acceptance

& \text{If}~u_{\alpha} <~\alpha, \

& \quad \text{Set}~q^k = q^*,~SS_{q^k} = SS_{q^*} \

& \text{Else} \

& \quad \text{Set}~q^k = q^{k-1},~SS_{q^k} = SS_{q^{k-1}}

Args:
Returns:
  • out_set (ParameterSet): If accept == 1, then latest DR set; Else, q^k=q^{k-1}

pymcmcstat.samplers.Metropolis module

Created on Thu Jan 18 10:30:29 2018

@author: prmiles

class pymcmcstat.samplers.Metropolis.Metropolis[source]

Pseudo-Algorithm:

  1. Sample z_k \sim N(0,1)
  2. Construct candidate q^* = q^{k-1} + Rz_k
  3. Compute
    \quad SS_{q^*} = \sum_{i=1}^N[v_i-f_i(q^*)]^2
  4. Compute
    \quad \alpha = \min\Big(1, e^{[SS_{q^*} - SS_{q^{k-1}}]/(2\sigma^2_{k-1})}\Big)
  5. If u_{\alpha} <~\alpha,
    Set q^k = q^*,~SS_{q^k} = SS_{q^*}
    Else
    Set q^k = q^{k-1},~SS_{q^k} = SS_{q^{k-1}}
Attributes:
classmethod evaluate_likelihood_function(ss1, ss2, sigma2, newprior, oldprior)[source]

Evaluate likelihood function:

\alpha = \exp\Big[-0.5\Big(\sum\Big(\frac{ SS_{q^*} - SS_{q^{k-1}} }{ \sigma_{k-1}^2 }\Big) + p_1 - p_2\Big)\Big]

Args:
  • ss1 (ndarray): SS error from proposed candidate, q^*
  • ss2 (ndarray): SS error from previous sample point, q^{k-1}
  • sigma2 (ndarray): Error variance estimate from previous sample point, \sigma_{k-1}^2
  • newprior (ndarray): Prior for proposal candidate
  • oldprior (ndarray): Prior for previous sample
Returns:
  • alpha (float): Result of likelihood function
run_metropolis_step(old_set, parameters, R, prior_object, sos_object)[source]

Run Metropolis step.

Args:
Returns:
  • accept (int): 0 - reject, 1 - accept
  • newset (ParameterSet): Features of q^*
  • outbound (int): 1 - rejected due to sampling outside of parameter bounds
  • npar_sample_from_normal (ndarray): Latet random sample points
classmethod unpack_set(parset)[source]

Unpack parameter set

Args:
Returns:
  • theta (ndarray): Value of sampled model parameters
  • ss (ndarray): Sum-of-squares error using sampled value
  • prior (ndarray): Value of prior
  • sigma2 (ndarray): Error variance

pymcmcstat.samplers.SamplingMethods module

Created on Thu Jan 18 10:11:40 2018

@author: prmiles

class pymcmcstat.samplers.SamplingMethods.SamplingMethods[source]

Metropolis sampling methods.

Attributes:

pymcmcstat.samplers.utilities module

Created on Wed Jun 27 12:07:11 2018

Utility functions used by different samplers

@author: prmiles

pymcmcstat.samplers.utilities.acceptance_test(alpha)[source]

Run standard acceptance test

& \text{If}~u_{\alpha} <~\alpha, \

& \quad \text{Set}~q^k = q^*,~SS_{q^k} = SS_{q^*} \

& \text{Else} \

& \quad \text{Set}~q^k = q^{k-1},~SS_{q^k} = SS_{q^{k-1}}

Args:
  • alpha (float): Result of likelihood function
Returns:
  • accept (bool): False - reject, True - accept
pymcmcstat.samplers.utilities.is_sample_outside_bounds(theta, lower_limits, upper_limits)[source]

Check whether proposal value is outside parameter limits

Args:
  • theta (ndarray): Value of sampled model parameters
  • lower_limits (ndarray): Lower limits
  • upper_limits (ndarray): Upper limits
Returns:
  • outsidebounds (bool): True -> Outside of parameter limits
pymcmcstat.samplers.utilities.sample_candidate_from_gaussian_proposal(npar, oldpar, R)[source]

Sample candidate from Gaussian proposal distribution

Args:
  • npar (int): Number of parameters being samples
  • oldpar (ndarray): q^{k-1} Old parameter set.
  • R (ndarray): Cholesky decomposition of parameter covariance matrix.
Returns:
  • newpar (ndarray): q^* - candidate
  • npar_sample_from_sample (ndarray): Sampled values from normal distibution: N(0,1).
pymcmcstat.samplers.utilities.set_outside_bounds(next_set)[source]

Assign set features based on being outside bounds

Args:
Returns: