Module mcmc
MCMC class for Bayes inference with emcee package using an Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler [1].
- Foreman-Mackey, D., Hogg, D.W., Lang, D. and Goodman, J., 2013. emcee: the MCMC hammer. Publications of the Astronomical Society of the Pacific, 125(925), p.306. https://emcee.readthedocs.io/en/stable/
Author: Farid Mohammadi, M.Sc. E-Mail: farid.mohammadi@iws.uni-stuttgart.de Department of Hydromechanics and Modelling of Hydrosystems (LH2) Institute for Modelling Hydraulic and Environmental Systems (IWS), University of Stuttgart, www.iws.uni-stuttgart.de/lh2/ Pfaffenwaldring 61 70569 Stuttgart
Created on Wed Jun 3 2020
Expand source code
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
MCMC class for Bayes inference with emcee package using an Affine Invariant
Markov chain Monte Carlo (MCMC) Ensemble sampler [1].
1. Foreman-Mackey, D., Hogg, D.W., Lang, D. and Goodman, J., 2013. emcee: the
MCMC hammer. Publications of the Astronomical Society of the Pacific, 125(925),
p.306. https://emcee.readthedocs.io/en/stable/
Author: Farid Mohammadi, M.Sc.
E-Mail: farid.mohammadi@iws.uni-stuttgart.de
Department of Hydromechanics and Modelling of Hydrosystems (LH2)
Institute for Modelling Hydraulic and Environmental Systems (IWS), University
of Stuttgart, www.iws.uni-stuttgart.de/lh2/
Pfaffenwaldring 61
70569 Stuttgart
Created on Wed Jun 3 2020
"""
import os
import numpy as np
import emcee
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
import multiprocessing
import scipy.stats as st
from scipy.linalg import cholesky as chol
import warnings
import shutil
os.environ["OMP_NUM_THREADS"] = "1"
class MCMC:
"""
A class for bayesian inference using a Markov-Chain Monte-Carlo Sampler.
"""
def __init__(self, BayesOpts):
self.BayesOpts = BayesOpts
def run_sampler(self, observation, total_sigma2):
BayesObj = self.BayesOpts
MetaModel = BayesObj.MetaModel
Model = MetaModel.ModelObj
Discrepancy = self.BayesOpts.Discrepancy
n_cpus = Model.n_cpus
priorDist = MetaModel.ExpDesign.JDist
ndim = MetaModel.n_params
self.counter = 0
output_dir = f'Outputs_Bayes_{Model.name}_{self.BayesOpts.name}'
if not os.path.exists(output_dir):
os.makedirs(output_dir)
self.observation = observation
self.total_sigma2 = total_sigma2
# Unpack mcmc parameters given to BayesObj.mcmc_params
self.initsamples = None
self.nwalkers = 100
self.nburn = 200
self.nsteps = 100000
self.moves = None
self.mp = False
self.verbose = False
# Extract initial samples
if 'init_samples' in BayesObj.mcmc_params:
self.initsamples = BayesObj.mcmc_params['init_samples']
if isinstance(self.initsamples, pd.DataFrame):
self.initsamples = self.initsamples.values
# Extract number of steps per walker
if 'n_steps' in BayesObj.mcmc_params:
self.nsteps = int(BayesObj.mcmc_params['n_steps'])
# Extract number of walkers (chains)
if 'n_walkers' in BayesObj.mcmc_params:
self.nwalkers = int(BayesObj.mcmc_params['n_walkers'])
# Extract moves
if 'moves' in BayesObj.mcmc_params:
self.moves = BayesObj.mcmc_params['moves']
# Extract multiprocessing
if 'multiprocessing' in BayesObj.mcmc_params:
self.mp = BayesObj.mcmc_params['multiprocessing']
# Extract verbose
if 'verbose' in BayesObj.mcmc_params:
self.verbose = BayesObj.mcmc_params['verbose']
# Set initial samples
np.random.seed(0)
if self.initsamples is None:
try:
initsamples = priorDist.sample(self.nwalkers).T
except:
# when aPCE selected - gaussian kernel distribution
inputSamples = MetaModel.ExpDesign.raw_data.T
random_indices = np.random.choice(
len(inputSamples), size=self.nwalkers, replace=False
)
initsamples = inputSamples[random_indices]
else:
if self.initsamples.ndim == 1:
# When MAL is given.
theta = self.initsamples
initsamples = [theta + 1e-1*np.multiply(
np.random.randn(ndim), theta) for i in
range(self.nwalkers)]
else:
# Pick samples based on a uniform dist between min and max of
# each dim
initsamples = np.zeros((self.nwalkers, ndim))
bound_tuples = []
for idx_dim in range(ndim):
lower = np.min(self.initsamples[:, idx_dim])
upper = np.max(self.initsamples[:, idx_dim])
bound_tuples.append((lower, upper))
dist = st.uniform(loc=lower, scale=upper-lower)
initsamples[:, idx_dim] = dist.rvs(size=self.nwalkers)
# Update lower and upper
MetaModel.ExpDesign.bound_tuples = bound_tuples
# Check if sigma^2 needs to be inferred
if Discrepancy.opt_sigma != 'B':
sigma2_samples = Discrepancy.get_sample(self.nwalkers)
# Update initsamples
initsamples = np.hstack((initsamples, sigma2_samples))
# Update ndim
ndim = initsamples.shape[1]
# Discrepancy bound
disc_bound_tuple = Discrepancy.ExpDesign.bound_tuples
# Update bound_tuples
MetaModel.ExpDesign.bound_tuples += disc_bound_tuple
print("\n>>>> Bayesian inference with MCMC for "
f"{self.BayesOpts.name} started. <<<<<<")
# Set up the backend
filename = f"{output_dir}/emcee_sampler.h5"
backend = emcee.backends.HDFBackend(filename)
# Clear the backend in case the file already exists
backend.reset(self.nwalkers, ndim)
# Define emcee sampler
# Here we'll set up the computation. emcee combines multiple "walkers",
# each of which is its own MCMC chain. The number of trace results will
# be nwalkers * nsteps.
if self.mp:
# Run in parallel
if n_cpus is None:
n_cpus = multiprocessing.cpu_count()
with multiprocessing.Pool(n_cpus) as pool:
sampler = emcee.EnsembleSampler(
self.nwalkers, ndim, self.log_posterior, moves=self.moves,
pool=pool, backend=backend
)
# Check if a burn-in phase is needed!
if self.initsamples is None:
# Burn-in
print("\n Burn-in period is starting:")
pos = sampler.run_mcmc(
initsamples, self.nburn, progress=True
)
# Reset sampler
sampler.reset()
pos = pos.coords
else:
pos = initsamples
# Production run
print("\n Production run is starting:")
pos, prob, state = sampler.run_mcmc(
pos, self.nsteps, progress=True
)
else:
# Run in series and monitor the convergence
sampler = emcee.EnsembleSampler(
self.nwalkers, ndim, self.log_posterior, moves=self.moves,
backend=backend, vectorize=True
)
# Check if a burn-in phase is needed!
if self.initsamples is None:
# Burn-in
print("\n Burn-in period is starting:")
pos = sampler.run_mcmc(
initsamples, self.nburn, progress=True
)
# Reset sampler
sampler.reset()
pos = pos.coords
else:
pos = initsamples
# Production run
print("\n Production run is starting:")
# Track how the average autocorrelation time estimate changes
autocorrIdx = 0
autocorr = np.empty(self.nsteps)
tauold = np.inf
autocorreverynsteps = 50
# sample step by step using the generator sampler.sample
for sample in sampler.sample(pos,
iterations=self.nsteps,
tune=True,
progress=True):
# only check convergence every autocorreverynsteps steps
if sampler.iteration % autocorreverynsteps:
continue
# Train model discrepancy/error
if hasattr(BayesObj, 'errorModel') and BayesObj.errorModel \
and not sampler.iteration % 3 * autocorreverynsteps:
try:
self.error_MetaModel = self.train_error_model(sampler)
except:
pass
# Print the current mean acceptance fraction
if self.verbose:
print("\nStep: {}".format(sampler.iteration))
acc_fr = np.mean(sampler.acceptance_fraction)
print(f"Mean acceptance fraction: {acc_fr:.3f}")
# compute the autocorrelation time so far
# using tol=0 means that we'll always get an estimate even if
# it isn't trustworthy
tau = sampler.get_autocorr_time(tol=0)
# average over walkers
autocorr[autocorrIdx] = np.nanmean(tau)
autocorrIdx += 1
# output current autocorrelation estimate
if self.verbose:
print(f"Mean autocorr time estimate: {np.nanmean(tau):.3f}")
list_gr = np.round(self.gelman_rubin(sampler.chain), 3)
print("Gelman-Rubin Test*: ", list_gr)
# check convergence
converged = np.all(tau*autocorreverynsteps < sampler.iteration)
converged &= np.all(np.abs(tauold - tau) / tau < 0.01)
converged &= np.all(self.gelman_rubin(sampler.chain) < 1.1)
if converged:
break
tauold = tau
# Posterior diagnostics
try:
tau = sampler.get_autocorr_time(tol=0)
except emcee.autocorr.AutocorrError:
tau = 5
if all(np.isnan(tau)):
tau = 5
burnin = int(2*np.nanmax(tau))
thin = int(0.5*np.nanmin(tau)) if int(0.5*np.nanmin(tau)) != 0 else 1
finalsamples = sampler.get_chain(discard=burnin, flat=True, thin=thin)
acc_fr = np.nanmean(sampler.acceptance_fraction)
list_gr = np.round(self.gelman_rubin(sampler.chain[:, burnin:]), 3)
# Print summary
print('\n')
print('-'*15 + 'Posterior diagnostics' + '-'*15)
print(f"mean auto-correlation time: {np.nanmean(tau):.3f}")
print(f"thin: {thin}")
print(f"burn-in: {burnin}")
print(f"flat chain shape: {finalsamples.shape}")
print(f"Mean acceptance fraction: {acc_fr:.3f}")
print("Gelman-Rubin Test*: ", list_gr)
print("\n* This value must lay between 0.234 and 0.5.")
print("* These values must be smaller than 1.1.")
print('-'*50)
print(f"\n>>>> Bayesian inference with MCMC for {self.BayesOpts.name} "
"successfully completed. <<<<<<\n")
# Extract parameter names and their prior ranges
par_names = MetaModel.ExpDesign.par_names
if Discrepancy.opt_sigma != 'B':
for i in range(len(Discrepancy.InputDisc.Marginals)):
par_names.append(Discrepancy.InputDisc.Marginals[i].name)
params_range = MetaModel.ExpDesign.bound_tuples
# Plot traces
if self.verbose and self.nsteps < 10000:
pdf = PdfPages(output_dir+'/traceplots.pdf')
fig = plt.figure()
for parIdx in range(ndim):
# Set up the axes with gridspec
fig = plt.figure()
grid = plt.GridSpec(4, 4, hspace=0.2, wspace=0.2)
main_ax = fig.add_subplot(grid[:-1, :3])
y_hist = fig.add_subplot(grid[:-1, -1], xticklabels=[],
sharey=main_ax)
for i in range(self.nwalkers):
samples = sampler.chain[i, :, parIdx]
main_ax.plot(samples, '-')
# histogram on the attached axes
y_hist.hist(samples[burnin:], 40, histtype='stepfilled',
orientation='horizontal', color='gray')
main_ax.set_ylim(params_range[parIdx])
main_ax.set_title('traceplot for ' + par_names[parIdx])
main_ax.set_xlabel('step number')
# save the current figure
pdf.savefig(fig, bbox_inches='tight')
# Destroy the current plot
plt.clf()
pdf.close()
# plot development of autocorrelation estimate
if not self.mp:
fig1 = plt.figure()
steps = autocorreverynsteps*np.arange(1, autocorrIdx+1)
taus = autocorr[:autocorrIdx]
plt.plot(steps, steps / autocorreverynsteps, "--k")
plt.plot(steps, taus)
plt.xlim(0, steps.max())
plt.ylim(0, np.nanmax(taus)+0.1*(np.nanmax(taus)-np.nanmin(taus)))
plt.xlabel("number of steps")
plt.ylabel(r"mean $\hat{\tau}$")
fig1.savefig(f"{output_dir}/autocorrelation_time.pdf",
bbox_inches='tight')
# logml_dict = self.marginal_llk_emcee(sampler, self.nburn, logp=None,
# maxiter=5000)
# print('\nThe Bridge Sampling Estimation is "
# f"{logml_dict['logml']:.5f}.')
# # Posterior-based expectation of posterior probablity
# postExpPostLikelihoods = np.mean(sampler.get_log_prob(flat=True)
# [self.nburn*self.nwalkers:])
# # Posterior-based expectation of prior densities
# postExpPrior = np.mean(self.log_prior(emcee_trace.T))
# # Posterior-based expectation of likelihoods
# postExpLikelihoods_emcee = postExpPostLikelihoods - postExpPrior
# # Calculate Kullback-Leibler Divergence
# KLD_emcee = postExpLikelihoods_emcee - logml_dict['logml']
# print("Kullback-Leibler divergence: %.5f"%KLD_emcee)
# # Information Entropy based on Entropy paper Eq. 38
# infEntropy_emcee = logml_dict['logml'] - postExpPrior -
# postExpLikelihoods_emcee
# print("Information Entropy: %.5f" %infEntropy_emcee)
Posterior_df = pd.DataFrame(finalsamples, columns=par_names)
return Posterior_df
# -------------------------------------------------------------------------
def log_prior(self, theta):
"""
Calculates the log prior likelihood for the given parameter set(s)
theta.
Parameters
----------
theta : array of shape (n_samples, n_params)
Parameter sets, i.e. proposals of MCMC chains.
Returns
-------
logprior: float or array of shape n_samples
Log prior likelihood. If theta has only one row, a single value is
returned otherwise an array.
"""
MetaModel = self.BayesOpts.MetaModel
Discrepancy = self.BayesOpts.Discrepancy
# Find the number of sigma2 parameters
if Discrepancy.opt_sigma != 'B':
disc_bound_tuples = Discrepancy.ExpDesign.bound_tuples
disc_marginals = Discrepancy.ExpDesign.InputObj.Marginals
disc_prior_space = Discrepancy.ExpDesign.prior_space
n_sigma2 = len(disc_bound_tuples)
else:
n_sigma2 = -len(theta)
prior_dist = MetaModel.ExpDesign.prior_space
params_range = MetaModel.ExpDesign.bound_tuples
theta = theta if theta.ndim != 1 else theta.reshape((1, -1))
nsamples = theta.shape[0]
logprior = -np.inf*np.ones(nsamples)
for i in range(nsamples):
# Check if the sample is within the parameters' range
if self._check_ranges(theta[i], params_range):
# Check if all dists are uniform, if yes priors are equal.
if all(MetaModel.input_obj.Marginals[i].dist_type == 'uniform'
for i in range(MetaModel.n_params)):
logprior[i] = 0.0
else:
logprior[i] = np.log(
prior_dist.pdf(theta[i, :-n_sigma2].T)
)
# Check if bias term needs to be inferred
if Discrepancy.opt_sigma != 'B':
if self._check_ranges(theta[i, -n_sigma2:],
disc_bound_tuples):
if all('unif' in disc_marginals[i].dist_type for i in
range(Discrepancy.ExpDesign.ndim)):
logprior[i] = 0.0
else:
logprior[i] += np.log(
disc_prior_space.pdf(theta[i, -n_sigma2:])
)
if nsamples == 1:
return logprior[0]
else:
return logprior
# -------------------------------------------------------------------------
def log_likelihood(self, theta):
"""
Computes likelihood p(y|theta, obs) of the performance of the
(meta-)model in reproducing the observation data.
Parameters
----------
theta : array of shape (n_samples, n_params)
Parameter set, i.e. proposals of the MCMC chains.
Returns
-------
log_like : array of shape (n_samples)
Log likelihood.
"""
BayesOpts = self.BayesOpts
MetaModel = BayesOpts.MetaModel
Discrepancy = self.BayesOpts.Discrepancy
# Find the number of sigma2 parameters
if Discrepancy.opt_sigma != 'B':
disc_bound_tuples = Discrepancy.ExpDesign.bound_tuples
n_sigma2 = len(disc_bound_tuples)
else:
n_sigma2 = -len(theta)
# Check if bias term needs to be inferred
if Discrepancy.opt_sigma != 'B':
sigma2 = theta[:, -n_sigma2:]
theta = theta[:, :-n_sigma2]
else:
sigma2 = None
theta = theta if theta.ndim != 1 else theta.reshape((1, -1))
# Evaluate Model/MetaModel at theta
mean_pred, BayesOpts._std_pce_prior_pred = self.eval_model(theta)
# Surrogate model's error using RMSE of test data
surrError = MetaModel.rmse if hasattr(MetaModel, 'rmse') else None
# Likelihood
log_like = BayesOpts.normpdf(
mean_pred, self.observation, self.total_sigma2, sigma2,
std=surrError
)
return log_like
# -------------------------------------------------------------------------
def log_posterior(self, theta):
"""
Computes the posterior likelihood p(theta| obs) for the given
parameterset.
Parameters
----------
theta : array of shape (n_samples, n_params)
Parameter set, i.e. proposals of the MCMC chains.
Returns
-------
log_like : array of shape (n_samples)
Log posterior likelihood.
"""
nsamples = 1 if theta.ndim == 1 else theta.shape[0]
if nsamples == 1:
if self.log_prior(theta) == -np.inf:
return -np.inf
else:
# Compute log prior
log_prior = self.log_prior(theta)
# Compute log Likelihood
log_likelihood = self.log_likelihood(theta)
return log_prior + log_likelihood
else:
# Compute log prior
log_prior = self.log_prior(theta)
# Initialize log_likelihood
log_likelihood = -np.inf*np.ones(nsamples)
# find the indices for -inf sets
non_inf_idx = np.where(log_prior != -np.inf)[0]
# Compute loLikelihoods
if non_inf_idx.size != 0:
log_likelihood[non_inf_idx] = self.log_likelihood(
theta[non_inf_idx]
)
return log_prior + log_likelihood
# -------------------------------------------------------------------------
def eval_model(self, theta):
"""
Evaluates the (meta-) model at the given theta.
Parameters
----------
theta : array of shape (n_samples, n_params)
Parameter set, i.e. proposals of the MCMC chains.
Returns
-------
mean_pred : dict
Mean model prediction.
std_pred : dict
Std of model prediction.
"""
BayesObj = self.BayesOpts
MetaModel = BayesObj.MetaModel
Model = MetaModel.ModelObj
if BayesObj.emulator:
# Evaluate the MetaModel
mean_pred, std_pred = MetaModel.eval_metamodel(samples=theta)
else:
# Evaluate the origModel
# Prepare the function
mean_pred, std_pred = dict(), dict()
output_names = Model.Output.names
try:
filename = Model.py_file
function = getattr(__import__(filename), filename)
# Run function for theta
model_outs = function(theta.reshape(1, -1))
except:
model_outs, _ = Model.run_model_parallel(
theta[0], prevRun_No=self.counter,
key_str='_MCMC', mp=False)
# Save outputs in respective dicts
for varIdx, var in enumerate(output_names):
mean_pred[var] = model_outs[var]
std_pred[var] = np.zeros((mean_pred[var].shape))
# Remove the folder
shutil.rmtree(f"{Model.Name}_MCMC_{self.counter+1}")
# Add one to the counter
self.counter += 1
if hasattr(self, 'error_MetaModel') and BayesObj.error_model:
meanPred, stdPred = self.error_MetaModel.eval_model_error(
BayesObj.BiasInputs, mean_pred
)
return mean_pred, std_pred
# -------------------------------------------------------------------------
def train_error_model(self, sampler):
"""
Trains an error model using a Gaussian Process Regression.
Parameters
----------
sampler : obj
emcee sampler.
Returns
-------
error_MetaModel : obj
A error model.
"""
BayesObj = self.BayesOpts
MetaModel = BayesObj.MetaModel
# Prepare the poster samples
try:
tau = sampler.get_autocorr_time(tol=0)
except emcee.autocorr.AutocorrError:
tau = 5
if all(np.isnan(tau)):
tau = 5
burnin = int(2*np.nanmax(tau))
thin = int(0.5*np.nanmin(tau)) if int(0.5*np.nanmin(tau)) != 0 else 1
finalsamples = sampler.get_chain(discard=burnin, flat=True, thin=thin)
posterior = finalsamples[:, :MetaModel.n_params]
# Select posterior mean as MAP
map_theta = posterior.mean(axis=0).reshape((1, MetaModel.n_params))
# MAP_theta = st.mode(Posterior_df,axis=0)[0]
# Evaluate the (meta-)model at the MAP
y_map, y_std_map = MetaModel.eval_metamodel(samples=map_theta)
# Train a GPR meta-model using MAP
error_MetaModel = MetaModel.create_model_error(
BayesObj.BiasInputs, y_map, name='Calib')
return error_MetaModel
# -------------------------------------------------------------------------
def gelman_rubin(self, chain, return_var=False):
"""
The potential scale reduction factor (PSRF) defined by the variance
within one chain, W, with the variance between chains B.
Both variances are combined in a weighted sum to obtain an estimate of
the variance of a parameter θ.The square root of the ratio of this
estimates variance to the within chain variance is called the potential
scale reduction.
For a well converged chain it should approach 1. Values greater than
typically 1.1 indicate that the chains have not yet fully converged.
Source: http://joergdietrich.github.io/emcee-convergence.html
https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py
Parameters
----------
chain : array (n_walkers, n_steps, n_params)
DESCRIPTION.
Returns
-------
R_hat : float
The Gelman-Robin values.
"""
m_chains, n_iters = chain.shape[:2]
# Calculate between-chain variance
θb = np.mean(chain, axis=1)
θbb = np.mean(θb, axis=0)
B_over_n = ((θbb - θb)**2).sum(axis=0)
B_over_n /= (m_chains - 1)
# Calculate within-chain variances
ssq = np.var(chain, axis=1, ddof=1)
W = np.mean(ssq, axis=0)
# (over) estimate of variance
var_θ = W * (n_iters - 1) / n_iters + B_over_n
if return_var:
return var_θ
else:
# The square root of the ratio of this estimates variance to the
# within chain variance
R_hat = np.sqrt(var_θ / W)
return R_hat
# -------------------------------------------------------------------------
def marginal_llk_emcee(self, sampler, nburn=None, logp=None, maxiter=1000):
"""
The Bridge Sampling Estimator of the Marginal Likelihood based on
https://gist.github.com/junpenglao/4d2669d69ddfe1d788318264cdcf0583
Parameters
----------
sampler : TYPE
MultiTrace, result of MCMC run.
nburn : int, optional
Number of burn-in step. The default is None.
logp : TYPE, optional
Model Log-probability function. The default is None.
maxiter : int, optional
Maximum number of iterations. The default is 1000.
Returns
-------
marg_llk : dict
Estimated Marginal log-Likelihood.
"""
r0, tol1, tol2 = 0.5, 1e-10, 1e-4
if logp is None:
logp = sampler.log_prob_fn
# Split the samples into two parts
# Use the first 50% for fiting the proposal distribution
# and the second 50% in the iterative scheme.
if nburn is None:
mtrace = sampler.chain
else:
mtrace = sampler.chain[:, nburn:, :]
nchain, len_trace, nrofVars = mtrace.shape
N1_ = len_trace // 2
N1 = N1_*nchain
N2 = len_trace*nchain - N1
samples_4_fit = np.zeros((nrofVars, N1))
samples_4_iter = np.zeros((nrofVars, N2))
effective_n = np.zeros((nrofVars))
# matrix with already transformed samples
for var in range(nrofVars):
# for fitting the proposal
x = mtrace[:, :N1_, var]
samples_4_fit[var, :] = x.flatten()
# for the iterative scheme
x2 = mtrace[:, N1_:, var]
samples_4_iter[var, :] = x2.flatten()
# effective sample size of samples_4_iter, scalar
effective_n[var] = self._my_ESS(x2)
# median effective sample size (scalar)
neff = np.median(effective_n)
# get mean & covariance matrix and generate samples from proposal
m = np.mean(samples_4_fit, axis=1)
V = np.cov(samples_4_fit)
L = chol(V, lower=True)
# Draw N2 samples from the proposal distribution
gen_samples = m[:, None] + np.dot(
L, st.norm.rvs(0, 1, size=samples_4_iter.shape)
)
# Evaluate proposal distribution for posterior & generated samples
q12 = st.multivariate_normal.logpdf(samples_4_iter.T, m, V)
q22 = st.multivariate_normal.logpdf(gen_samples.T, m, V)
# Evaluate unnormalized posterior for posterior & generated samples
q11 = logp(samples_4_iter.T)
q21 = logp(gen_samples.T)
# Run iterative scheme:
tmp = self._iterative_scheme(
N1, N2, q11, q12, q21, q22, r0, neff, tol1, maxiter, 'r'
)
if ~np.isfinite(tmp['logml']):
warnings.warn(
"logml could not be estimated within maxiter, rerunning with "
"adjusted starting value. Estimate might be more variable than"
" usual.")
# use geometric mean as starting value
r0_2 = np.sqrt(tmp['r_vals'][-2]*tmp['r_vals'][-1])
tmp = self._iterative_scheme(
q11, q12, q21, q22, r0_2, neff, tol2, maxiter, 'logml'
)
marg_llk = dict(
logml=tmp['logml'], niter=tmp['niter'], method="normal",
q11=q11, q12=q12, q21=q21, q22=q22
)
return marg_llk
# -------------------------------------------------------------------------
def _iterative_scheme(self, N1, N2, q11, q12, q21, q22, r0, neff, tol,
maxiter, criterion):
"""
Iterative scheme as proposed in Meng and Wong (1996) to estimate the
marginal likelihood
Parameters
----------
N1 : TYPE
DESCRIPTION.
N2 : TYPE
DESCRIPTION.
q11 : TYPE
DESCRIPTION.
q12 : TYPE
DESCRIPTION.
q21 : TYPE
DESCRIPTION.
q22 : TYPE
DESCRIPTION.
r0 : TYPE
DESCRIPTION.
neff : TYPE
DESCRIPTION.
tol : TYPE
DESCRIPTION.
maxiter : TYPE
DESCRIPTION.
criterion : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
l1 = q11 - q12
l2 = q21 - q22
# To increase numerical stability,
# subtracting the median of l1 from l1 & l2 later
lstar = np.median(l1)
s1 = neff/(neff + N2)
s2 = N2/(neff + N2)
r = r0
r_vals = [r]
logml = np.log(r) + lstar
criterion_val = 1 + tol
i = 0
while (i <= maxiter) & (criterion_val > tol):
rold = r
logmlold = logml
numi = np.exp(l2 - lstar)/(s1 * np.exp(l2 - lstar) + s2 * r)
deni = 1/(s1 * np.exp(l1 - lstar) + s2 * r)
if np.sum(~np.isfinite(numi))+np.sum(~np.isfinite(deni)) > 0:
warnings.warn(
"""Infinite value in iterative scheme, returning NaN.
Try rerunning with more samples.""")
r = (N1/N2) * np.sum(numi)/np.sum(deni)
r_vals.append(r)
logml = np.log(r) + lstar
i += 1
if criterion == 'r':
criterion_val = np.abs((r - rold)/r)
elif criterion == 'logml':
criterion_val = np.abs((logml - logmlold)/logml)
if i >= maxiter:
return dict(logml=np.NaN, niter=i, r_vals=np.asarray(r_vals))
else:
return dict(logml=logml, niter=i)
# -------------------------------------------------------------------------
def _my_ESS(self, x):
"""
Compute the effective sample size of estimand of interest.
Vectorised implementation.
https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py
Parameters
----------
x : array of shape (n_walkers, n_steps)
MCMC Samples.
Returns
-------
int
Effective sample size.
"""
m_chains, n_iters = x.shape
def variogram(t):
variogram = ((x[:, t:] - x[:, :(n_iters - t)])**2).sum()
variogram /= (m_chains * (n_iters - t))
return variogram
post_var = self.gelman_rubin(x, return_var=True)
t = 1
rho = np.ones(n_iters)
negative_autocorr = False
# Iterate until the sum of consecutive estimates of autocorrelation is
# negative
while not negative_autocorr and (t < n_iters):
rho[t] = 1 - variogram(t) / (2 * post_var)
if not t % 2:
negative_autocorr = sum(rho[t-1:t+1]) < 0
t += 1
return int(m_chains*n_iters / (1 + 2*rho[1:t].sum()))
# -------------------------------------------------------------------------
def _check_ranges(self, theta, ranges):
"""
This function checks if theta lies in the given ranges.
Parameters
----------
theta : array
Proposed parameter set.
ranges : nested list
List of the praremeter ranges.
Returns
-------
c : bool
If it lies in the given range, it return True else False.
"""
c = True
# traverse in the list1
for i, bounds in enumerate(ranges):
x = theta[i]
# condition check
if x < bounds[0] or x > bounds[1]:
c = False
return c
return cClasses
class MCMC (BayesOpts)-
A class for bayesian inference using a Markov-Chain Monte-Carlo Sampler.
Expand source code
class MCMC: """ A class for bayesian inference using a Markov-Chain Monte-Carlo Sampler. """ def __init__(self, BayesOpts): self.BayesOpts = BayesOpts def run_sampler(self, observation, total_sigma2): BayesObj = self.BayesOpts MetaModel = BayesObj.MetaModel Model = MetaModel.ModelObj Discrepancy = self.BayesOpts.Discrepancy n_cpus = Model.n_cpus priorDist = MetaModel.ExpDesign.JDist ndim = MetaModel.n_params self.counter = 0 output_dir = f'Outputs_Bayes_{Model.name}_{self.BayesOpts.name}' if not os.path.exists(output_dir): os.makedirs(output_dir) self.observation = observation self.total_sigma2 = total_sigma2 # Unpack mcmc parameters given to BayesObj.mcmc_params self.initsamples = None self.nwalkers = 100 self.nburn = 200 self.nsteps = 100000 self.moves = None self.mp = False self.verbose = False # Extract initial samples if 'init_samples' in BayesObj.mcmc_params: self.initsamples = BayesObj.mcmc_params['init_samples'] if isinstance(self.initsamples, pd.DataFrame): self.initsamples = self.initsamples.values # Extract number of steps per walker if 'n_steps' in BayesObj.mcmc_params: self.nsteps = int(BayesObj.mcmc_params['n_steps']) # Extract number of walkers (chains) if 'n_walkers' in BayesObj.mcmc_params: self.nwalkers = int(BayesObj.mcmc_params['n_walkers']) # Extract moves if 'moves' in BayesObj.mcmc_params: self.moves = BayesObj.mcmc_params['moves'] # Extract multiprocessing if 'multiprocessing' in BayesObj.mcmc_params: self.mp = BayesObj.mcmc_params['multiprocessing'] # Extract verbose if 'verbose' in BayesObj.mcmc_params: self.verbose = BayesObj.mcmc_params['verbose'] # Set initial samples np.random.seed(0) if self.initsamples is None: try: initsamples = priorDist.sample(self.nwalkers).T except: # when aPCE selected - gaussian kernel distribution inputSamples = MetaModel.ExpDesign.raw_data.T random_indices = np.random.choice( len(inputSamples), size=self.nwalkers, replace=False ) initsamples = inputSamples[random_indices] else: if self.initsamples.ndim == 1: # When MAL is given. theta = self.initsamples initsamples = [theta + 1e-1*np.multiply( np.random.randn(ndim), theta) for i in range(self.nwalkers)] else: # Pick samples based on a uniform dist between min and max of # each dim initsamples = np.zeros((self.nwalkers, ndim)) bound_tuples = [] for idx_dim in range(ndim): lower = np.min(self.initsamples[:, idx_dim]) upper = np.max(self.initsamples[:, idx_dim]) bound_tuples.append((lower, upper)) dist = st.uniform(loc=lower, scale=upper-lower) initsamples[:, idx_dim] = dist.rvs(size=self.nwalkers) # Update lower and upper MetaModel.ExpDesign.bound_tuples = bound_tuples # Check if sigma^2 needs to be inferred if Discrepancy.opt_sigma != 'B': sigma2_samples = Discrepancy.get_sample(self.nwalkers) # Update initsamples initsamples = np.hstack((initsamples, sigma2_samples)) # Update ndim ndim = initsamples.shape[1] # Discrepancy bound disc_bound_tuple = Discrepancy.ExpDesign.bound_tuples # Update bound_tuples MetaModel.ExpDesign.bound_tuples += disc_bound_tuple print("\n>>>> Bayesian inference with MCMC for " f"{self.BayesOpts.name} started. <<<<<<") # Set up the backend filename = f"{output_dir}/emcee_sampler.h5" backend = emcee.backends.HDFBackend(filename) # Clear the backend in case the file already exists backend.reset(self.nwalkers, ndim) # Define emcee sampler # Here we'll set up the computation. emcee combines multiple "walkers", # each of which is its own MCMC chain. The number of trace results will # be nwalkers * nsteps. if self.mp: # Run in parallel if n_cpus is None: n_cpus = multiprocessing.cpu_count() with multiprocessing.Pool(n_cpus) as pool: sampler = emcee.EnsembleSampler( self.nwalkers, ndim, self.log_posterior, moves=self.moves, pool=pool, backend=backend ) # Check if a burn-in phase is needed! if self.initsamples is None: # Burn-in print("\n Burn-in period is starting:") pos = sampler.run_mcmc( initsamples, self.nburn, progress=True ) # Reset sampler sampler.reset() pos = pos.coords else: pos = initsamples # Production run print("\n Production run is starting:") pos, prob, state = sampler.run_mcmc( pos, self.nsteps, progress=True ) else: # Run in series and monitor the convergence sampler = emcee.EnsembleSampler( self.nwalkers, ndim, self.log_posterior, moves=self.moves, backend=backend, vectorize=True ) # Check if a burn-in phase is needed! if self.initsamples is None: # Burn-in print("\n Burn-in period is starting:") pos = sampler.run_mcmc( initsamples, self.nburn, progress=True ) # Reset sampler sampler.reset() pos = pos.coords else: pos = initsamples # Production run print("\n Production run is starting:") # Track how the average autocorrelation time estimate changes autocorrIdx = 0 autocorr = np.empty(self.nsteps) tauold = np.inf autocorreverynsteps = 50 # sample step by step using the generator sampler.sample for sample in sampler.sample(pos, iterations=self.nsteps, tune=True, progress=True): # only check convergence every autocorreverynsteps steps if sampler.iteration % autocorreverynsteps: continue # Train model discrepancy/error if hasattr(BayesObj, 'errorModel') and BayesObj.errorModel \ and not sampler.iteration % 3 * autocorreverynsteps: try: self.error_MetaModel = self.train_error_model(sampler) except: pass # Print the current mean acceptance fraction if self.verbose: print("\nStep: {}".format(sampler.iteration)) acc_fr = np.mean(sampler.acceptance_fraction) print(f"Mean acceptance fraction: {acc_fr:.3f}") # compute the autocorrelation time so far # using tol=0 means that we'll always get an estimate even if # it isn't trustworthy tau = sampler.get_autocorr_time(tol=0) # average over walkers autocorr[autocorrIdx] = np.nanmean(tau) autocorrIdx += 1 # output current autocorrelation estimate if self.verbose: print(f"Mean autocorr time estimate: {np.nanmean(tau):.3f}") list_gr = np.round(self.gelman_rubin(sampler.chain), 3) print("Gelman-Rubin Test*: ", list_gr) # check convergence converged = np.all(tau*autocorreverynsteps < sampler.iteration) converged &= np.all(np.abs(tauold - tau) / tau < 0.01) converged &= np.all(self.gelman_rubin(sampler.chain) < 1.1) if converged: break tauold = tau # Posterior diagnostics try: tau = sampler.get_autocorr_time(tol=0) except emcee.autocorr.AutocorrError: tau = 5 if all(np.isnan(tau)): tau = 5 burnin = int(2*np.nanmax(tau)) thin = int(0.5*np.nanmin(tau)) if int(0.5*np.nanmin(tau)) != 0 else 1 finalsamples = sampler.get_chain(discard=burnin, flat=True, thin=thin) acc_fr = np.nanmean(sampler.acceptance_fraction) list_gr = np.round(self.gelman_rubin(sampler.chain[:, burnin:]), 3) # Print summary print('\n') print('-'*15 + 'Posterior diagnostics' + '-'*15) print(f"mean auto-correlation time: {np.nanmean(tau):.3f}") print(f"thin: {thin}") print(f"burn-in: {burnin}") print(f"flat chain shape: {finalsamples.shape}") print(f"Mean acceptance fraction: {acc_fr:.3f}") print("Gelman-Rubin Test*: ", list_gr) print("\n* This value must lay between 0.234 and 0.5.") print("* These values must be smaller than 1.1.") print('-'*50) print(f"\n>>>> Bayesian inference with MCMC for {self.BayesOpts.name} " "successfully completed. <<<<<<\n") # Extract parameter names and their prior ranges par_names = MetaModel.ExpDesign.par_names if Discrepancy.opt_sigma != 'B': for i in range(len(Discrepancy.InputDisc.Marginals)): par_names.append(Discrepancy.InputDisc.Marginals[i].name) params_range = MetaModel.ExpDesign.bound_tuples # Plot traces if self.verbose and self.nsteps < 10000: pdf = PdfPages(output_dir+'/traceplots.pdf') fig = plt.figure() for parIdx in range(ndim): # Set up the axes with gridspec fig = plt.figure() grid = plt.GridSpec(4, 4, hspace=0.2, wspace=0.2) main_ax = fig.add_subplot(grid[:-1, :3]) y_hist = fig.add_subplot(grid[:-1, -1], xticklabels=[], sharey=main_ax) for i in range(self.nwalkers): samples = sampler.chain[i, :, parIdx] main_ax.plot(samples, '-') # histogram on the attached axes y_hist.hist(samples[burnin:], 40, histtype='stepfilled', orientation='horizontal', color='gray') main_ax.set_ylim(params_range[parIdx]) main_ax.set_title('traceplot for ' + par_names[parIdx]) main_ax.set_xlabel('step number') # save the current figure pdf.savefig(fig, bbox_inches='tight') # Destroy the current plot plt.clf() pdf.close() # plot development of autocorrelation estimate if not self.mp: fig1 = plt.figure() steps = autocorreverynsteps*np.arange(1, autocorrIdx+1) taus = autocorr[:autocorrIdx] plt.plot(steps, steps / autocorreverynsteps, "--k") plt.plot(steps, taus) plt.xlim(0, steps.max()) plt.ylim(0, np.nanmax(taus)+0.1*(np.nanmax(taus)-np.nanmin(taus))) plt.xlabel("number of steps") plt.ylabel(r"mean $\hat{\tau}$") fig1.savefig(f"{output_dir}/autocorrelation_time.pdf", bbox_inches='tight') # logml_dict = self.marginal_llk_emcee(sampler, self.nburn, logp=None, # maxiter=5000) # print('\nThe Bridge Sampling Estimation is " # f"{logml_dict['logml']:.5f}.') # # Posterior-based expectation of posterior probablity # postExpPostLikelihoods = np.mean(sampler.get_log_prob(flat=True) # [self.nburn*self.nwalkers:]) # # Posterior-based expectation of prior densities # postExpPrior = np.mean(self.log_prior(emcee_trace.T)) # # Posterior-based expectation of likelihoods # postExpLikelihoods_emcee = postExpPostLikelihoods - postExpPrior # # Calculate Kullback-Leibler Divergence # KLD_emcee = postExpLikelihoods_emcee - logml_dict['logml'] # print("Kullback-Leibler divergence: %.5f"%KLD_emcee) # # Information Entropy based on Entropy paper Eq. 38 # infEntropy_emcee = logml_dict['logml'] - postExpPrior - # postExpLikelihoods_emcee # print("Information Entropy: %.5f" %infEntropy_emcee) Posterior_df = pd.DataFrame(finalsamples, columns=par_names) return Posterior_df # ------------------------------------------------------------------------- def log_prior(self, theta): """ Calculates the log prior likelihood for the given parameter set(s) theta. Parameters ---------- theta : array of shape (n_samples, n_params) Parameter sets, i.e. proposals of MCMC chains. Returns ------- logprior: float or array of shape n_samples Log prior likelihood. If theta has only one row, a single value is returned otherwise an array. """ MetaModel = self.BayesOpts.MetaModel Discrepancy = self.BayesOpts.Discrepancy # Find the number of sigma2 parameters if Discrepancy.opt_sigma != 'B': disc_bound_tuples = Discrepancy.ExpDesign.bound_tuples disc_marginals = Discrepancy.ExpDesign.InputObj.Marginals disc_prior_space = Discrepancy.ExpDesign.prior_space n_sigma2 = len(disc_bound_tuples) else: n_sigma2 = -len(theta) prior_dist = MetaModel.ExpDesign.prior_space params_range = MetaModel.ExpDesign.bound_tuples theta = theta if theta.ndim != 1 else theta.reshape((1, -1)) nsamples = theta.shape[0] logprior = -np.inf*np.ones(nsamples) for i in range(nsamples): # Check if the sample is within the parameters' range if self._check_ranges(theta[i], params_range): # Check if all dists are uniform, if yes priors are equal. if all(MetaModel.input_obj.Marginals[i].dist_type == 'uniform' for i in range(MetaModel.n_params)): logprior[i] = 0.0 else: logprior[i] = np.log( prior_dist.pdf(theta[i, :-n_sigma2].T) ) # Check if bias term needs to be inferred if Discrepancy.opt_sigma != 'B': if self._check_ranges(theta[i, -n_sigma2:], disc_bound_tuples): if all('unif' in disc_marginals[i].dist_type for i in range(Discrepancy.ExpDesign.ndim)): logprior[i] = 0.0 else: logprior[i] += np.log( disc_prior_space.pdf(theta[i, -n_sigma2:]) ) if nsamples == 1: return logprior[0] else: return logprior # ------------------------------------------------------------------------- def log_likelihood(self, theta): """ Computes likelihood p(y|theta, obs) of the performance of the (meta-)model in reproducing the observation data. Parameters ---------- theta : array of shape (n_samples, n_params) Parameter set, i.e. proposals of the MCMC chains. Returns ------- log_like : array of shape (n_samples) Log likelihood. """ BayesOpts = self.BayesOpts MetaModel = BayesOpts.MetaModel Discrepancy = self.BayesOpts.Discrepancy # Find the number of sigma2 parameters if Discrepancy.opt_sigma != 'B': disc_bound_tuples = Discrepancy.ExpDesign.bound_tuples n_sigma2 = len(disc_bound_tuples) else: n_sigma2 = -len(theta) # Check if bias term needs to be inferred if Discrepancy.opt_sigma != 'B': sigma2 = theta[:, -n_sigma2:] theta = theta[:, :-n_sigma2] else: sigma2 = None theta = theta if theta.ndim != 1 else theta.reshape((1, -1)) # Evaluate Model/MetaModel at theta mean_pred, BayesOpts._std_pce_prior_pred = self.eval_model(theta) # Surrogate model's error using RMSE of test data surrError = MetaModel.rmse if hasattr(MetaModel, 'rmse') else None # Likelihood log_like = BayesOpts.normpdf( mean_pred, self.observation, self.total_sigma2, sigma2, std=surrError ) return log_like # ------------------------------------------------------------------------- def log_posterior(self, theta): """ Computes the posterior likelihood p(theta| obs) for the given parameterset. Parameters ---------- theta : array of shape (n_samples, n_params) Parameter set, i.e. proposals of the MCMC chains. Returns ------- log_like : array of shape (n_samples) Log posterior likelihood. """ nsamples = 1 if theta.ndim == 1 else theta.shape[0] if nsamples == 1: if self.log_prior(theta) == -np.inf: return -np.inf else: # Compute log prior log_prior = self.log_prior(theta) # Compute log Likelihood log_likelihood = self.log_likelihood(theta) return log_prior + log_likelihood else: # Compute log prior log_prior = self.log_prior(theta) # Initialize log_likelihood log_likelihood = -np.inf*np.ones(nsamples) # find the indices for -inf sets non_inf_idx = np.where(log_prior != -np.inf)[0] # Compute loLikelihoods if non_inf_idx.size != 0: log_likelihood[non_inf_idx] = self.log_likelihood( theta[non_inf_idx] ) return log_prior + log_likelihood # ------------------------------------------------------------------------- def eval_model(self, theta): """ Evaluates the (meta-) model at the given theta. Parameters ---------- theta : array of shape (n_samples, n_params) Parameter set, i.e. proposals of the MCMC chains. Returns ------- mean_pred : dict Mean model prediction. std_pred : dict Std of model prediction. """ BayesObj = self.BayesOpts MetaModel = BayesObj.MetaModel Model = MetaModel.ModelObj if BayesObj.emulator: # Evaluate the MetaModel mean_pred, std_pred = MetaModel.eval_metamodel(samples=theta) else: # Evaluate the origModel # Prepare the function mean_pred, std_pred = dict(), dict() output_names = Model.Output.names try: filename = Model.py_file function = getattr(__import__(filename), filename) # Run function for theta model_outs = function(theta.reshape(1, -1)) except: model_outs, _ = Model.run_model_parallel( theta[0], prevRun_No=self.counter, key_str='_MCMC', mp=False) # Save outputs in respective dicts for varIdx, var in enumerate(output_names): mean_pred[var] = model_outs[var] std_pred[var] = np.zeros((mean_pred[var].shape)) # Remove the folder shutil.rmtree(f"{Model.Name}_MCMC_{self.counter+1}") # Add one to the counter self.counter += 1 if hasattr(self, 'error_MetaModel') and BayesObj.error_model: meanPred, stdPred = self.error_MetaModel.eval_model_error( BayesObj.BiasInputs, mean_pred ) return mean_pred, std_pred # ------------------------------------------------------------------------- def train_error_model(self, sampler): """ Trains an error model using a Gaussian Process Regression. Parameters ---------- sampler : obj emcee sampler. Returns ------- error_MetaModel : obj A error model. """ BayesObj = self.BayesOpts MetaModel = BayesObj.MetaModel # Prepare the poster samples try: tau = sampler.get_autocorr_time(tol=0) except emcee.autocorr.AutocorrError: tau = 5 if all(np.isnan(tau)): tau = 5 burnin = int(2*np.nanmax(tau)) thin = int(0.5*np.nanmin(tau)) if int(0.5*np.nanmin(tau)) != 0 else 1 finalsamples = sampler.get_chain(discard=burnin, flat=True, thin=thin) posterior = finalsamples[:, :MetaModel.n_params] # Select posterior mean as MAP map_theta = posterior.mean(axis=0).reshape((1, MetaModel.n_params)) # MAP_theta = st.mode(Posterior_df,axis=0)[0] # Evaluate the (meta-)model at the MAP y_map, y_std_map = MetaModel.eval_metamodel(samples=map_theta) # Train a GPR meta-model using MAP error_MetaModel = MetaModel.create_model_error( BayesObj.BiasInputs, y_map, name='Calib') return error_MetaModel # ------------------------------------------------------------------------- def gelman_rubin(self, chain, return_var=False): """ The potential scale reduction factor (PSRF) defined by the variance within one chain, W, with the variance between chains B. Both variances are combined in a weighted sum to obtain an estimate of the variance of a parameter θ.The square root of the ratio of this estimates variance to the within chain variance is called the potential scale reduction. For a well converged chain it should approach 1. Values greater than typically 1.1 indicate that the chains have not yet fully converged. Source: http://joergdietrich.github.io/emcee-convergence.html https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py Parameters ---------- chain : array (n_walkers, n_steps, n_params) DESCRIPTION. Returns ------- R_hat : float The Gelman-Robin values. """ m_chains, n_iters = chain.shape[:2] # Calculate between-chain variance θb = np.mean(chain, axis=1) θbb = np.mean(θb, axis=0) B_over_n = ((θbb - θb)**2).sum(axis=0) B_over_n /= (m_chains - 1) # Calculate within-chain variances ssq = np.var(chain, axis=1, ddof=1) W = np.mean(ssq, axis=0) # (over) estimate of variance var_θ = W * (n_iters - 1) / n_iters + B_over_n if return_var: return var_θ else: # The square root of the ratio of this estimates variance to the # within chain variance R_hat = np.sqrt(var_θ / W) return R_hat # ------------------------------------------------------------------------- def marginal_llk_emcee(self, sampler, nburn=None, logp=None, maxiter=1000): """ The Bridge Sampling Estimator of the Marginal Likelihood based on https://gist.github.com/junpenglao/4d2669d69ddfe1d788318264cdcf0583 Parameters ---------- sampler : TYPE MultiTrace, result of MCMC run. nburn : int, optional Number of burn-in step. The default is None. logp : TYPE, optional Model Log-probability function. The default is None. maxiter : int, optional Maximum number of iterations. The default is 1000. Returns ------- marg_llk : dict Estimated Marginal log-Likelihood. """ r0, tol1, tol2 = 0.5, 1e-10, 1e-4 if logp is None: logp = sampler.log_prob_fn # Split the samples into two parts # Use the first 50% for fiting the proposal distribution # and the second 50% in the iterative scheme. if nburn is None: mtrace = sampler.chain else: mtrace = sampler.chain[:, nburn:, :] nchain, len_trace, nrofVars = mtrace.shape N1_ = len_trace // 2 N1 = N1_*nchain N2 = len_trace*nchain - N1 samples_4_fit = np.zeros((nrofVars, N1)) samples_4_iter = np.zeros((nrofVars, N2)) effective_n = np.zeros((nrofVars)) # matrix with already transformed samples for var in range(nrofVars): # for fitting the proposal x = mtrace[:, :N1_, var] samples_4_fit[var, :] = x.flatten() # for the iterative scheme x2 = mtrace[:, N1_:, var] samples_4_iter[var, :] = x2.flatten() # effective sample size of samples_4_iter, scalar effective_n[var] = self._my_ESS(x2) # median effective sample size (scalar) neff = np.median(effective_n) # get mean & covariance matrix and generate samples from proposal m = np.mean(samples_4_fit, axis=1) V = np.cov(samples_4_fit) L = chol(V, lower=True) # Draw N2 samples from the proposal distribution gen_samples = m[:, None] + np.dot( L, st.norm.rvs(0, 1, size=samples_4_iter.shape) ) # Evaluate proposal distribution for posterior & generated samples q12 = st.multivariate_normal.logpdf(samples_4_iter.T, m, V) q22 = st.multivariate_normal.logpdf(gen_samples.T, m, V) # Evaluate unnormalized posterior for posterior & generated samples q11 = logp(samples_4_iter.T) q21 = logp(gen_samples.T) # Run iterative scheme: tmp = self._iterative_scheme( N1, N2, q11, q12, q21, q22, r0, neff, tol1, maxiter, 'r' ) if ~np.isfinite(tmp['logml']): warnings.warn( "logml could not be estimated within maxiter, rerunning with " "adjusted starting value. Estimate might be more variable than" " usual.") # use geometric mean as starting value r0_2 = np.sqrt(tmp['r_vals'][-2]*tmp['r_vals'][-1]) tmp = self._iterative_scheme( q11, q12, q21, q22, r0_2, neff, tol2, maxiter, 'logml' ) marg_llk = dict( logml=tmp['logml'], niter=tmp['niter'], method="normal", q11=q11, q12=q12, q21=q21, q22=q22 ) return marg_llk # ------------------------------------------------------------------------- def _iterative_scheme(self, N1, N2, q11, q12, q21, q22, r0, neff, tol, maxiter, criterion): """ Iterative scheme as proposed in Meng and Wong (1996) to estimate the marginal likelihood Parameters ---------- N1 : TYPE DESCRIPTION. N2 : TYPE DESCRIPTION. q11 : TYPE DESCRIPTION. q12 : TYPE DESCRIPTION. q21 : TYPE DESCRIPTION. q22 : TYPE DESCRIPTION. r0 : TYPE DESCRIPTION. neff : TYPE DESCRIPTION. tol : TYPE DESCRIPTION. maxiter : TYPE DESCRIPTION. criterion : TYPE DESCRIPTION. Returns ------- TYPE DESCRIPTION. """ l1 = q11 - q12 l2 = q21 - q22 # To increase numerical stability, # subtracting the median of l1 from l1 & l2 later lstar = np.median(l1) s1 = neff/(neff + N2) s2 = N2/(neff + N2) r = r0 r_vals = [r] logml = np.log(r) + lstar criterion_val = 1 + tol i = 0 while (i <= maxiter) & (criterion_val > tol): rold = r logmlold = logml numi = np.exp(l2 - lstar)/(s1 * np.exp(l2 - lstar) + s2 * r) deni = 1/(s1 * np.exp(l1 - lstar) + s2 * r) if np.sum(~np.isfinite(numi))+np.sum(~np.isfinite(deni)) > 0: warnings.warn( """Infinite value in iterative scheme, returning NaN. Try rerunning with more samples.""") r = (N1/N2) * np.sum(numi)/np.sum(deni) r_vals.append(r) logml = np.log(r) + lstar i += 1 if criterion == 'r': criterion_val = np.abs((r - rold)/r) elif criterion == 'logml': criterion_val = np.abs((logml - logmlold)/logml) if i >= maxiter: return dict(logml=np.NaN, niter=i, r_vals=np.asarray(r_vals)) else: return dict(logml=logml, niter=i) # ------------------------------------------------------------------------- def _my_ESS(self, x): """ Compute the effective sample size of estimand of interest. Vectorised implementation. https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py Parameters ---------- x : array of shape (n_walkers, n_steps) MCMC Samples. Returns ------- int Effective sample size. """ m_chains, n_iters = x.shape def variogram(t): variogram = ((x[:, t:] - x[:, :(n_iters - t)])**2).sum() variogram /= (m_chains * (n_iters - t)) return variogram post_var = self.gelman_rubin(x, return_var=True) t = 1 rho = np.ones(n_iters) negative_autocorr = False # Iterate until the sum of consecutive estimates of autocorrelation is # negative while not negative_autocorr and (t < n_iters): rho[t] = 1 - variogram(t) / (2 * post_var) if not t % 2: negative_autocorr = sum(rho[t-1:t+1]) < 0 t += 1 return int(m_chains*n_iters / (1 + 2*rho[1:t].sum())) # ------------------------------------------------------------------------- def _check_ranges(self, theta, ranges): """ This function checks if theta lies in the given ranges. Parameters ---------- theta : array Proposed parameter set. ranges : nested list List of the praremeter ranges. Returns ------- c : bool If it lies in the given range, it return True else False. """ c = True # traverse in the list1 for i, bounds in enumerate(ranges): x = theta[i] # condition check if x < bounds[0] or x > bounds[1]: c = False return c return cMethods
def eval_model(self, theta)-
Evaluates the (meta-) model at the given theta.
Parameters
theta:arrayofshape (n_samples, n_params)- Parameter set, i.e. proposals of the MCMC chains.
Returns
mean_pred:dict- Mean model prediction.
std_pred:dict- Std of model prediction.
Expand source code
def eval_model(self, theta): """ Evaluates the (meta-) model at the given theta. Parameters ---------- theta : array of shape (n_samples, n_params) Parameter set, i.e. proposals of the MCMC chains. Returns ------- mean_pred : dict Mean model prediction. std_pred : dict Std of model prediction. """ BayesObj = self.BayesOpts MetaModel = BayesObj.MetaModel Model = MetaModel.ModelObj if BayesObj.emulator: # Evaluate the MetaModel mean_pred, std_pred = MetaModel.eval_metamodel(samples=theta) else: # Evaluate the origModel # Prepare the function mean_pred, std_pred = dict(), dict() output_names = Model.Output.names try: filename = Model.py_file function = getattr(__import__(filename), filename) # Run function for theta model_outs = function(theta.reshape(1, -1)) except: model_outs, _ = Model.run_model_parallel( theta[0], prevRun_No=self.counter, key_str='_MCMC', mp=False) # Save outputs in respective dicts for varIdx, var in enumerate(output_names): mean_pred[var] = model_outs[var] std_pred[var] = np.zeros((mean_pred[var].shape)) # Remove the folder shutil.rmtree(f"{Model.Name}_MCMC_{self.counter+1}") # Add one to the counter self.counter += 1 if hasattr(self, 'error_MetaModel') and BayesObj.error_model: meanPred, stdPred = self.error_MetaModel.eval_model_error( BayesObj.BiasInputs, mean_pred ) return mean_pred, std_pred def gelman_rubin(self, chain, return_var=False)-
The potential scale reduction factor (PSRF) defined by the variance within one chain, W, with the variance between chains B. Both variances are combined in a weighted sum to obtain an estimate of the variance of a parameter θ.The square root of the ratio of this estimates variance to the within chain variance is called the potential scale reduction. For a well converged chain it should approach 1. Values greater than typically 1.1 indicate that the chains have not yet fully converged.
Source: http://joergdietrich.github.io/emcee-convergence.html
https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py
Parameters
chain:array (n_walkers, n_steps, n_params)- DESCRIPTION.
Returns
R_hat:float- The Gelman-Robin values.
Expand source code
def gelman_rubin(self, chain, return_var=False): """ The potential scale reduction factor (PSRF) defined by the variance within one chain, W, with the variance between chains B. Both variances are combined in a weighted sum to obtain an estimate of the variance of a parameter θ.The square root of the ratio of this estimates variance to the within chain variance is called the potential scale reduction. For a well converged chain it should approach 1. Values greater than typically 1.1 indicate that the chains have not yet fully converged. Source: http://joergdietrich.github.io/emcee-convergence.html https://github.com/jwalton3141/jwalton3141.github.io/blob/master/assets/posts/ESS/rwmh.py Parameters ---------- chain : array (n_walkers, n_steps, n_params) DESCRIPTION. Returns ------- R_hat : float The Gelman-Robin values. """ m_chains, n_iters = chain.shape[:2] # Calculate between-chain variance θb = np.mean(chain, axis=1) θbb = np.mean(θb, axis=0) B_over_n = ((θbb - θb)**2).sum(axis=0) B_over_n /= (m_chains - 1) # Calculate within-chain variances ssq = np.var(chain, axis=1, ddof=1) W = np.mean(ssq, axis=0) # (over) estimate of variance var_θ = W * (n_iters - 1) / n_iters + B_over_n if return_var: return var_θ else: # The square root of the ratio of this estimates variance to the # within chain variance R_hat = np.sqrt(var_θ / W) return R_hat def log_likelihood(self, theta)-
Computes likelihood p(y|theta, obs) of the performance of the (meta-)model in reproducing the observation data.
Parameters
theta:arrayofshape (n_samples, n_params)- Parameter set, i.e. proposals of the MCMC chains.
Returns
log_like:arrayofshape (n_samples)- Log likelihood.
Expand source code
def log_likelihood(self, theta): """ Computes likelihood p(y|theta, obs) of the performance of the (meta-)model in reproducing the observation data. Parameters ---------- theta : array of shape (n_samples, n_params) Parameter set, i.e. proposals of the MCMC chains. Returns ------- log_like : array of shape (n_samples) Log likelihood. """ BayesOpts = self.BayesOpts MetaModel = BayesOpts.MetaModel Discrepancy = self.BayesOpts.Discrepancy # Find the number of sigma2 parameters if Discrepancy.opt_sigma != 'B': disc_bound_tuples = Discrepancy.ExpDesign.bound_tuples n_sigma2 = len(disc_bound_tuples) else: n_sigma2 = -len(theta) # Check if bias term needs to be inferred if Discrepancy.opt_sigma != 'B': sigma2 = theta[:, -n_sigma2:] theta = theta[:, :-n_sigma2] else: sigma2 = None theta = theta if theta.ndim != 1 else theta.reshape((1, -1)) # Evaluate Model/MetaModel at theta mean_pred, BayesOpts._std_pce_prior_pred = self.eval_model(theta) # Surrogate model's error using RMSE of test data surrError = MetaModel.rmse if hasattr(MetaModel, 'rmse') else None # Likelihood log_like = BayesOpts.normpdf( mean_pred, self.observation, self.total_sigma2, sigma2, std=surrError ) return log_like def log_posterior(self, theta)-
Computes the posterior likelihood p(theta| obs) for the given parameterset.
Parameters
theta:arrayofshape (n_samples, n_params)- Parameter set, i.e. proposals of the MCMC chains.
Returns
log_like:arrayofshape (n_samples)- Log posterior likelihood.
Expand source code
def log_posterior(self, theta): """ Computes the posterior likelihood p(theta| obs) for the given parameterset. Parameters ---------- theta : array of shape (n_samples, n_params) Parameter set, i.e. proposals of the MCMC chains. Returns ------- log_like : array of shape (n_samples) Log posterior likelihood. """ nsamples = 1 if theta.ndim == 1 else theta.shape[0] if nsamples == 1: if self.log_prior(theta) == -np.inf: return -np.inf else: # Compute log prior log_prior = self.log_prior(theta) # Compute log Likelihood log_likelihood = self.log_likelihood(theta) return log_prior + log_likelihood else: # Compute log prior log_prior = self.log_prior(theta) # Initialize log_likelihood log_likelihood = -np.inf*np.ones(nsamples) # find the indices for -inf sets non_inf_idx = np.where(log_prior != -np.inf)[0] # Compute loLikelihoods if non_inf_idx.size != 0: log_likelihood[non_inf_idx] = self.log_likelihood( theta[non_inf_idx] ) return log_prior + log_likelihood def log_prior(self, theta)-
Calculates the log prior likelihood for the given parameter set(s) theta.
Parameters
theta:arrayofshape (n_samples, n_params)- Parameter sets, i.e. proposals of MCMC chains.
Returns
logprior:floatorarrayofshape n_samples- Log prior likelihood. If theta has only one row, a single value is returned otherwise an array.
Expand source code
def log_prior(self, theta): """ Calculates the log prior likelihood for the given parameter set(s) theta. Parameters ---------- theta : array of shape (n_samples, n_params) Parameter sets, i.e. proposals of MCMC chains. Returns ------- logprior: float or array of shape n_samples Log prior likelihood. If theta has only one row, a single value is returned otherwise an array. """ MetaModel = self.BayesOpts.MetaModel Discrepancy = self.BayesOpts.Discrepancy # Find the number of sigma2 parameters if Discrepancy.opt_sigma != 'B': disc_bound_tuples = Discrepancy.ExpDesign.bound_tuples disc_marginals = Discrepancy.ExpDesign.InputObj.Marginals disc_prior_space = Discrepancy.ExpDesign.prior_space n_sigma2 = len(disc_bound_tuples) else: n_sigma2 = -len(theta) prior_dist = MetaModel.ExpDesign.prior_space params_range = MetaModel.ExpDesign.bound_tuples theta = theta if theta.ndim != 1 else theta.reshape((1, -1)) nsamples = theta.shape[0] logprior = -np.inf*np.ones(nsamples) for i in range(nsamples): # Check if the sample is within the parameters' range if self._check_ranges(theta[i], params_range): # Check if all dists are uniform, if yes priors are equal. if all(MetaModel.input_obj.Marginals[i].dist_type == 'uniform' for i in range(MetaModel.n_params)): logprior[i] = 0.0 else: logprior[i] = np.log( prior_dist.pdf(theta[i, :-n_sigma2].T) ) # Check if bias term needs to be inferred if Discrepancy.opt_sigma != 'B': if self._check_ranges(theta[i, -n_sigma2:], disc_bound_tuples): if all('unif' in disc_marginals[i].dist_type for i in range(Discrepancy.ExpDesign.ndim)): logprior[i] = 0.0 else: logprior[i] += np.log( disc_prior_space.pdf(theta[i, -n_sigma2:]) ) if nsamples == 1: return logprior[0] else: return logprior def marginal_llk_emcee(self, sampler, nburn=None, logp=None, maxiter=1000)-
The Bridge Sampling Estimator of the Marginal Likelihood based on https://gist.github.com/junpenglao/4d2669d69ddfe1d788318264cdcf0583
Parameters
sampler:TYPE- MultiTrace, result of MCMC run.
nburn:int, optional- Number of burn-in step. The default is None.
logp:TYPE, optional- Model Log-probability function. The default is None.
maxiter:int, optional- Maximum number of iterations. The default is 1000.
Returns
marg_llk:dict- Estimated Marginal log-Likelihood.
Expand source code
def marginal_llk_emcee(self, sampler, nburn=None, logp=None, maxiter=1000): """ The Bridge Sampling Estimator of the Marginal Likelihood based on https://gist.github.com/junpenglao/4d2669d69ddfe1d788318264cdcf0583 Parameters ---------- sampler : TYPE MultiTrace, result of MCMC run. nburn : int, optional Number of burn-in step. The default is None. logp : TYPE, optional Model Log-probability function. The default is None. maxiter : int, optional Maximum number of iterations. The default is 1000. Returns ------- marg_llk : dict Estimated Marginal log-Likelihood. """ r0, tol1, tol2 = 0.5, 1e-10, 1e-4 if logp is None: logp = sampler.log_prob_fn # Split the samples into two parts # Use the first 50% for fiting the proposal distribution # and the second 50% in the iterative scheme. if nburn is None: mtrace = sampler.chain else: mtrace = sampler.chain[:, nburn:, :] nchain, len_trace, nrofVars = mtrace.shape N1_ = len_trace // 2 N1 = N1_*nchain N2 = len_trace*nchain - N1 samples_4_fit = np.zeros((nrofVars, N1)) samples_4_iter = np.zeros((nrofVars, N2)) effective_n = np.zeros((nrofVars)) # matrix with already transformed samples for var in range(nrofVars): # for fitting the proposal x = mtrace[:, :N1_, var] samples_4_fit[var, :] = x.flatten() # for the iterative scheme x2 = mtrace[:, N1_:, var] samples_4_iter[var, :] = x2.flatten() # effective sample size of samples_4_iter, scalar effective_n[var] = self._my_ESS(x2) # median effective sample size (scalar) neff = np.median(effective_n) # get mean & covariance matrix and generate samples from proposal m = np.mean(samples_4_fit, axis=1) V = np.cov(samples_4_fit) L = chol(V, lower=True) # Draw N2 samples from the proposal distribution gen_samples = m[:, None] + np.dot( L, st.norm.rvs(0, 1, size=samples_4_iter.shape) ) # Evaluate proposal distribution for posterior & generated samples q12 = st.multivariate_normal.logpdf(samples_4_iter.T, m, V) q22 = st.multivariate_normal.logpdf(gen_samples.T, m, V) # Evaluate unnormalized posterior for posterior & generated samples q11 = logp(samples_4_iter.T) q21 = logp(gen_samples.T) # Run iterative scheme: tmp = self._iterative_scheme( N1, N2, q11, q12, q21, q22, r0, neff, tol1, maxiter, 'r' ) if ~np.isfinite(tmp['logml']): warnings.warn( "logml could not be estimated within maxiter, rerunning with " "adjusted starting value. Estimate might be more variable than" " usual.") # use geometric mean as starting value r0_2 = np.sqrt(tmp['r_vals'][-2]*tmp['r_vals'][-1]) tmp = self._iterative_scheme( q11, q12, q21, q22, r0_2, neff, tol2, maxiter, 'logml' ) marg_llk = dict( logml=tmp['logml'], niter=tmp['niter'], method="normal", q11=q11, q12=q12, q21=q21, q22=q22 ) return marg_llk def run_sampler(self, observation, total_sigma2)-
Expand source code
def run_sampler(self, observation, total_sigma2): BayesObj = self.BayesOpts MetaModel = BayesObj.MetaModel Model = MetaModel.ModelObj Discrepancy = self.BayesOpts.Discrepancy n_cpus = Model.n_cpus priorDist = MetaModel.ExpDesign.JDist ndim = MetaModel.n_params self.counter = 0 output_dir = f'Outputs_Bayes_{Model.name}_{self.BayesOpts.name}' if not os.path.exists(output_dir): os.makedirs(output_dir) self.observation = observation self.total_sigma2 = total_sigma2 # Unpack mcmc parameters given to BayesObj.mcmc_params self.initsamples = None self.nwalkers = 100 self.nburn = 200 self.nsteps = 100000 self.moves = None self.mp = False self.verbose = False # Extract initial samples if 'init_samples' in BayesObj.mcmc_params: self.initsamples = BayesObj.mcmc_params['init_samples'] if isinstance(self.initsamples, pd.DataFrame): self.initsamples = self.initsamples.values # Extract number of steps per walker if 'n_steps' in BayesObj.mcmc_params: self.nsteps = int(BayesObj.mcmc_params['n_steps']) # Extract number of walkers (chains) if 'n_walkers' in BayesObj.mcmc_params: self.nwalkers = int(BayesObj.mcmc_params['n_walkers']) # Extract moves if 'moves' in BayesObj.mcmc_params: self.moves = BayesObj.mcmc_params['moves'] # Extract multiprocessing if 'multiprocessing' in BayesObj.mcmc_params: self.mp = BayesObj.mcmc_params['multiprocessing'] # Extract verbose if 'verbose' in BayesObj.mcmc_params: self.verbose = BayesObj.mcmc_params['verbose'] # Set initial samples np.random.seed(0) if self.initsamples is None: try: initsamples = priorDist.sample(self.nwalkers).T except: # when aPCE selected - gaussian kernel distribution inputSamples = MetaModel.ExpDesign.raw_data.T random_indices = np.random.choice( len(inputSamples), size=self.nwalkers, replace=False ) initsamples = inputSamples[random_indices] else: if self.initsamples.ndim == 1: # When MAL is given. theta = self.initsamples initsamples = [theta + 1e-1*np.multiply( np.random.randn(ndim), theta) for i in range(self.nwalkers)] else: # Pick samples based on a uniform dist between min and max of # each dim initsamples = np.zeros((self.nwalkers, ndim)) bound_tuples = [] for idx_dim in range(ndim): lower = np.min(self.initsamples[:, idx_dim]) upper = np.max(self.initsamples[:, idx_dim]) bound_tuples.append((lower, upper)) dist = st.uniform(loc=lower, scale=upper-lower) initsamples[:, idx_dim] = dist.rvs(size=self.nwalkers) # Update lower and upper MetaModel.ExpDesign.bound_tuples = bound_tuples # Check if sigma^2 needs to be inferred if Discrepancy.opt_sigma != 'B': sigma2_samples = Discrepancy.get_sample(self.nwalkers) # Update initsamples initsamples = np.hstack((initsamples, sigma2_samples)) # Update ndim ndim = initsamples.shape[1] # Discrepancy bound disc_bound_tuple = Discrepancy.ExpDesign.bound_tuples # Update bound_tuples MetaModel.ExpDesign.bound_tuples += disc_bound_tuple print("\n>>>> Bayesian inference with MCMC for " f"{self.BayesOpts.name} started. <<<<<<") # Set up the backend filename = f"{output_dir}/emcee_sampler.h5" backend = emcee.backends.HDFBackend(filename) # Clear the backend in case the file already exists backend.reset(self.nwalkers, ndim) # Define emcee sampler # Here we'll set up the computation. emcee combines multiple "walkers", # each of which is its own MCMC chain. The number of trace results will # be nwalkers * nsteps. if self.mp: # Run in parallel if n_cpus is None: n_cpus = multiprocessing.cpu_count() with multiprocessing.Pool(n_cpus) as pool: sampler = emcee.EnsembleSampler( self.nwalkers, ndim, self.log_posterior, moves=self.moves, pool=pool, backend=backend ) # Check if a burn-in phase is needed! if self.initsamples is None: # Burn-in print("\n Burn-in period is starting:") pos = sampler.run_mcmc( initsamples, self.nburn, progress=True ) # Reset sampler sampler.reset() pos = pos.coords else: pos = initsamples # Production run print("\n Production run is starting:") pos, prob, state = sampler.run_mcmc( pos, self.nsteps, progress=True ) else: # Run in series and monitor the convergence sampler = emcee.EnsembleSampler( self.nwalkers, ndim, self.log_posterior, moves=self.moves, backend=backend, vectorize=True ) # Check if a burn-in phase is needed! if self.initsamples is None: # Burn-in print("\n Burn-in period is starting:") pos = sampler.run_mcmc( initsamples, self.nburn, progress=True ) # Reset sampler sampler.reset() pos = pos.coords else: pos = initsamples # Production run print("\n Production run is starting:") # Track how the average autocorrelation time estimate changes autocorrIdx = 0 autocorr = np.empty(self.nsteps) tauold = np.inf autocorreverynsteps = 50 # sample step by step using the generator sampler.sample for sample in sampler.sample(pos, iterations=self.nsteps, tune=True, progress=True): # only check convergence every autocorreverynsteps steps if sampler.iteration % autocorreverynsteps: continue # Train model discrepancy/error if hasattr(BayesObj, 'errorModel') and BayesObj.errorModel \ and not sampler.iteration % 3 * autocorreverynsteps: try: self.error_MetaModel = self.train_error_model(sampler) except: pass # Print the current mean acceptance fraction if self.verbose: print("\nStep: {}".format(sampler.iteration)) acc_fr = np.mean(sampler.acceptance_fraction) print(f"Mean acceptance fraction: {acc_fr:.3f}") # compute the autocorrelation time so far # using tol=0 means that we'll always get an estimate even if # it isn't trustworthy tau = sampler.get_autocorr_time(tol=0) # average over walkers autocorr[autocorrIdx] = np.nanmean(tau) autocorrIdx += 1 # output current autocorrelation estimate if self.verbose: print(f"Mean autocorr time estimate: {np.nanmean(tau):.3f}") list_gr = np.round(self.gelman_rubin(sampler.chain), 3) print("Gelman-Rubin Test*: ", list_gr) # check convergence converged = np.all(tau*autocorreverynsteps < sampler.iteration) converged &= np.all(np.abs(tauold - tau) / tau < 0.01) converged &= np.all(self.gelman_rubin(sampler.chain) < 1.1) if converged: break tauold = tau # Posterior diagnostics try: tau = sampler.get_autocorr_time(tol=0) except emcee.autocorr.AutocorrError: tau = 5 if all(np.isnan(tau)): tau = 5 burnin = int(2*np.nanmax(tau)) thin = int(0.5*np.nanmin(tau)) if int(0.5*np.nanmin(tau)) != 0 else 1 finalsamples = sampler.get_chain(discard=burnin, flat=True, thin=thin) acc_fr = np.nanmean(sampler.acceptance_fraction) list_gr = np.round(self.gelman_rubin(sampler.chain[:, burnin:]), 3) # Print summary print('\n') print('-'*15 + 'Posterior diagnostics' + '-'*15) print(f"mean auto-correlation time: {np.nanmean(tau):.3f}") print(f"thin: {thin}") print(f"burn-in: {burnin}") print(f"flat chain shape: {finalsamples.shape}") print(f"Mean acceptance fraction: {acc_fr:.3f}") print("Gelman-Rubin Test*: ", list_gr) print("\n* This value must lay between 0.234 and 0.5.") print("* These values must be smaller than 1.1.") print('-'*50) print(f"\n>>>> Bayesian inference with MCMC for {self.BayesOpts.name} " "successfully completed. <<<<<<\n") # Extract parameter names and their prior ranges par_names = MetaModel.ExpDesign.par_names if Discrepancy.opt_sigma != 'B': for i in range(len(Discrepancy.InputDisc.Marginals)): par_names.append(Discrepancy.InputDisc.Marginals[i].name) params_range = MetaModel.ExpDesign.bound_tuples # Plot traces if self.verbose and self.nsteps < 10000: pdf = PdfPages(output_dir+'/traceplots.pdf') fig = plt.figure() for parIdx in range(ndim): # Set up the axes with gridspec fig = plt.figure() grid = plt.GridSpec(4, 4, hspace=0.2, wspace=0.2) main_ax = fig.add_subplot(grid[:-1, :3]) y_hist = fig.add_subplot(grid[:-1, -1], xticklabels=[], sharey=main_ax) for i in range(self.nwalkers): samples = sampler.chain[i, :, parIdx] main_ax.plot(samples, '-') # histogram on the attached axes y_hist.hist(samples[burnin:], 40, histtype='stepfilled', orientation='horizontal', color='gray') main_ax.set_ylim(params_range[parIdx]) main_ax.set_title('traceplot for ' + par_names[parIdx]) main_ax.set_xlabel('step number') # save the current figure pdf.savefig(fig, bbox_inches='tight') # Destroy the current plot plt.clf() pdf.close() # plot development of autocorrelation estimate if not self.mp: fig1 = plt.figure() steps = autocorreverynsteps*np.arange(1, autocorrIdx+1) taus = autocorr[:autocorrIdx] plt.plot(steps, steps / autocorreverynsteps, "--k") plt.plot(steps, taus) plt.xlim(0, steps.max()) plt.ylim(0, np.nanmax(taus)+0.1*(np.nanmax(taus)-np.nanmin(taus))) plt.xlabel("number of steps") plt.ylabel(r"mean $\hat{\tau}$") fig1.savefig(f"{output_dir}/autocorrelation_time.pdf", bbox_inches='tight') # logml_dict = self.marginal_llk_emcee(sampler, self.nburn, logp=None, # maxiter=5000) # print('\nThe Bridge Sampling Estimation is " # f"{logml_dict['logml']:.5f}.') # # Posterior-based expectation of posterior probablity # postExpPostLikelihoods = np.mean(sampler.get_log_prob(flat=True) # [self.nburn*self.nwalkers:]) # # Posterior-based expectation of prior densities # postExpPrior = np.mean(self.log_prior(emcee_trace.T)) # # Posterior-based expectation of likelihoods # postExpLikelihoods_emcee = postExpPostLikelihoods - postExpPrior # # Calculate Kullback-Leibler Divergence # KLD_emcee = postExpLikelihoods_emcee - logml_dict['logml'] # print("Kullback-Leibler divergence: %.5f"%KLD_emcee) # # Information Entropy based on Entropy paper Eq. 38 # infEntropy_emcee = logml_dict['logml'] - postExpPrior - # postExpLikelihoods_emcee # print("Information Entropy: %.5f" %infEntropy_emcee) Posterior_df = pd.DataFrame(finalsamples, columns=par_names) return Posterior_df def train_error_model(self, sampler)-
Trains an error model using a Gaussian Process Regression.
Parameters
sampler:obj- emcee sampler.
Returns
error_MetaModel:obj- A error model.
Expand source code
def train_error_model(self, sampler): """ Trains an error model using a Gaussian Process Regression. Parameters ---------- sampler : obj emcee sampler. Returns ------- error_MetaModel : obj A error model. """ BayesObj = self.BayesOpts MetaModel = BayesObj.MetaModel # Prepare the poster samples try: tau = sampler.get_autocorr_time(tol=0) except emcee.autocorr.AutocorrError: tau = 5 if all(np.isnan(tau)): tau = 5 burnin = int(2*np.nanmax(tau)) thin = int(0.5*np.nanmin(tau)) if int(0.5*np.nanmin(tau)) != 0 else 1 finalsamples = sampler.get_chain(discard=burnin, flat=True, thin=thin) posterior = finalsamples[:, :MetaModel.n_params] # Select posterior mean as MAP map_theta = posterior.mean(axis=0).reshape((1, MetaModel.n_params)) # MAP_theta = st.mode(Posterior_df,axis=0)[0] # Evaluate the (meta-)model at the MAP y_map, y_std_map = MetaModel.eval_metamodel(samples=map_theta) # Train a GPR meta-model using MAP error_MetaModel = MetaModel.create_model_error( BayesObj.BiasInputs, y_map, name='Calib') return error_MetaModel