pymcmcstat.samplers package¶

pymcmcstat.samplers.Adaptation module¶

Created on Thu Jan 18 11:14:11 2018

@author: prmiles

class pymcmcstat.samplers.Adaptation.Adaptation[source]¶

Bases: object

Adaptive Metropolis (AM) algorithm based on [HST+01].

Attributes:

run_adaptation()

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:

covariance (CovarianceProcedures): Updated covariance object

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:

x (ndarray): Chain segment
w (ndarray): Weights

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:

x (ndarray): Covariance matrix

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:

covariance (CovarianceProcedures): Covariance methods and variables

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:

options (SimulationOptions): Options for MCMC simulation

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]¶

Bases: object

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

Attributes:

run_delayed_rejection()
initialize_next_metropolis_step()
alphafun()

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, custom=None)[source]¶

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

Args:

old_set (ParameterSet): Features of $q^{k-1}$
new_set (ParameterSet): Features of $q^*$
RDR (ndarray): Cholesky decomposition of parameter covariance matrix for DR steps
ntry (int): Number of DR steps to perform until rejection
parameters (ModelParameters): Model parameters
invR (ndarray): Inverse Cholesky decomposition matrix
sosobj (SumOfSquares): Sum-of-Squares function
priorobj (PriorFunction): Prior function

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:

x1 (ParameterSet): Old set - $q^{k-1}$
x2 (ParameterSet): New set - $q^*$

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:

accept (int): 0 - reject, 1 - accept
old_set (ParameterSet): Features of $q^{k-1}$
next_set (ParameterSet): New proposal set

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]¶

Bases: object

Pseudo-Algorithm:

Sample $z_k \sim N(0,1)$

Construct candidate $q^* = q^{k-1} + Rz_k$

Compute
$\quad SS_{q^*} = \sum_{i=1}^N[v_i-f_i(q^*)]^2$

Compute
$\quad \alpha = \min\Big(1, e^{[SS_{q^*} - SS_{q^{k-1}}]/(2\sigma^2_{k-1})}\Big)$

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:

acceptance_test()
run_metropolis_step()
unpack_set()

classmethod calculate_posterior_ratio(ss1, ss2, sigma2, newprior, oldprior)[source]¶

Calculate acceptance ratio

$\alpha = \min\Big[1, \frac{\mathcal{L}(\nu_{obs}|q^*, \sigma_{k-1}^2)\pi_0(q^*)}{\mathcal{L}(\nu_{obs}|q^{k-1}, \sigma_{k-1}^2)\pi_0(q^{k-1})}\Big]$

where the Gaussian likelihood function is

$\mathcal{L}(\nu_{obs}|q, \sigma) = \exp\Big(-\frac{SS_q}{2\sigma}\Big)$

and Gaussian prior function is

$\pi_0(q) = \exp \Big[-\frac{1}{2}\Big(\frac{q - \mu_0}{\sigma_0}\Big)^2\Big].$

For the Gaussian likelihood and prior, this yields the acceptance ratio

$\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].$

For more details regarding the prior function, please refer to the PriorFunction class.

Note

The default behavior of the package is to use Gaussian likelihood and prior functions (as of v1.8.0). Future releases will expand the functionality to allow for alternative likelihood and prior definitions.

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, custom=None)[source]¶

Run Metropolis step.

Args:

old_set (ParameterSet): Features of $q^{k-1}$
parameters (ModelParameters): Model parameters
R (ndarray): Cholesky decomposition of parameter covariance matrix
priorobj (PriorFunction): Prior function
sosobj (SumOfSquares): Sum-of-Squares function

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:

parset (ParameterSet): Parameter set to unpack

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]¶

Bases: object

Metropolis sampling methods.

Attributes:

Metropolis
DelayedRejection
Adaptation

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.calculate_log_posterior_ratio(loglikestar, loglike, logpriorstar, logprior)[source]¶

Calculate log posterior ratio:

$\log(\alpha) = \min\Big[0, \log(\mathcal{L}(\nu_{obs}|q^*)) + \log(\pi_0(q^*)) - \log(\mathcal{L}(\nu_{obs}|q^{k-1})) - \log(\pi_0(q^{k-1}))\Big]$

For more details regarding the prior and likelihood functions, please refer to the PriorFunction and LikelihoodFunction class, respectively.

Note

The default behavior of the package is to use Gaussian likelihood and prior functions (as of v1.8.0). Future releases will expand the functionality to allow for alternative likelihood and prior definitions.

Args:

likestar (float): Likelihood from proposed candidate, $q^*$
like (float): Likelihood from previous sample point, $q^{k-1}$
priorstar (float): Prior for proposal candidate
prior (float): Prior for previous sample

Returns:

alpha (float): Acceptance ratio

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.log_posterior_ratio_acceptance_test(alpha)[source]¶

Run log posterior ratio acceptance test

Args:

alpha (float): Log posterior ratio

Returns:

accept (bool): False - reject, True - accept

pymcmcstat.samplers.utilities.posterior_ratio_acceptance_test(alpha)[source]¶

Run posterior ratio acceptance test

Args:

alpha (float): Posterior ratio

Returns:

accept (bool): False - reject, True - accept

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:

next_set (ParameterSet): $q^*$

Returns:

next_set (ParameterSet): $q^*$ with updated features
outbound (bool): True