Source code for snowdrop.src.numeric.bayes.mcmc

"""
Created on Fri Apr  5 13:09:34 2019

Curently implemented for linear models only.

@author: A.Goumilevski
"""

import os,sys
import numpy as np
import multiprocessing
from time import time
#from numeric.bayes.logp import logp
from snowdrop.src.misc.termcolor import cprint
from snowdrop.src.utils.distributions import pdf
from snowdrop.src.utils.distributions import getHyperParameters
from snowdrop.src.numeric.solver import linear_solver
from snowdrop.src.numeric.filters.filters import DK_Filter as dk_filter
from snowdrop.src.numeric.filters.filters import Durbin_Koopman_Non_Diffuse_Filter as dk_non_diffuse_filter 
from snowdrop.src.model.settings import SamplingAlgorithm
from snowdrop.src.model.settings import PriorAssumption,FilterAlgorithm
import snowdrop.src.numeric.solver.linear_solver as ls
from snowdrop.src.utils.prettyTable import PrettyTable

num_cores = multiprocessing.cpu_count()
num_chains = 3
PARALLEL = False # if True runs parallel sampling

count=0; nt=-1; n_=-1; n_shocks_=-1; ind_non_missing_=[]
model_=None; params=None; param_names=None; y0_=None; param_index=None; meas_shocks_=[]
obs_=None; Z_=None; Hm_=None; Qm_=None; pm_model=None; Pstar=None
pm_names=None; pm_params=None; priors=None; data_=None; est_shocks_names=[]; shocks_=[]
LARGE_NUMBER = 1.e10

path = os.path.dirname(os.path.abspath(__file__))
path_to_dir = os.path.join(path,'../../../graphs')

[docs] def prior_logp(p): """Function defining log of prior probability.""" global model_,params,param_names,param_index,count log_likelihood = 0 try: b = True for i,index in enumerate(param_index): name = param_names[index] prior = model_.priors[name] distr = prior['distribution'] pars = np.copy(prior['parameters']) x = p[i] lower_bound = float(pars[1]) upper_bound = float(pars[2]) if x < lower_bound or x > upper_bound: return -LARGE_NUMBER - count # Get distribution hyperparameters if not distr.endswith("_hp"): pars[3],pars[4] = getHyperParameters(distr=distr,mean=pars[3],std=pars[4]) likelihood,b = pdf(distr,x,pars) if not b: return -LARGE_NUMBER - count if np.isnan(likelihood) or np.isinf(likelihood): log_likelihood = -LARGE_NUMBER - count elif likelihood > 0: log_likelihood += np.log(max(1.e-10,likelihood)) # Find standard deviations of shocks npar = len(param_index) for i,name in enumerate(est_shocks_names): par = p[i+npar] prior = model_.priors[name] distr = prior['distribution'] pars = prior['parameters'] lower_bound = float(pars[1]) upper_bound = float(pars[2]) if par < lower_bound or par > upper_bound: return -LARGE_NUMBER - count else: # Get distribution hyperparameters if not distr.endswith("_hp"): pars[3],pars[4] = getHyperParameters(distr=distr,mean=pars[3],std=pars[4]) x,b = pdf(distr,par,pars) if not b: return -LARGE_NUMBER - count log_likelihood += np.log(max(1.e-10,x)) except: log_likelihood = -LARGE_NUMBER - count return log_likelihood
[docs] def likelihood_logp(parameters,stds): """Function defining log of likelihood.""" global model_,count,n_,nt,n_shocks_,Pstar,shocks_,meas_shocks_ global y0_,Qm_,obs_,Z_,Hm_,ind_non_missing_ from snowdrop.src.numeric.solver.util import getCovarianceMatrix # Expected value of outcome log_likelihood = 0 K = None; S = None; Sinv = None calib = {} for i,name in enumerate(est_shocks_names): calib[name] = stds[i] # Set values of covariance matrix Qm,Hm = getCovarianceMatrix(Qm_,Hm_,calib,shocks_,meas_shocks_) # Solve linear model model_.solved=False linear_solver.solve(model=model_,p=parameters,steady_state=np.zeros(n_),suppress_warnings=True) # State transition matrix F = model_.linear_model["A"] F1 = F[:n_,:n_] # Array of constants C = model_.linear_model["C"] C = C[:n_] # Matrix of coefficients of shocks R = model_.linear_model["R"] R = R[:n_] Q = 0 for i in range(1+model_.max_lead_shock-model_.min_lag_shock): R1 = R[:,i*n_shocks_:(1+i)*n_shocks_] Q += np.real(R1 @ Qm_ @ R1.T) bUnivariate = True y = yt = np.copy(y0_) # Get model likelihood by applying Kalman filter for t in range (nt): # Initialize covariance matrices if t == 0: P = np.copy(Pstar) # Apply Kalman filter to find log-likelyhood if not model_.FILTER is None and model_.FILTER.value == FilterAlgorithm.Non_Diffuse_Filter.value: yt,y,residual,P,_,K,S,Sinv,bUnivariate,loglk = \ dk_non_diffuse_filter(x=yt,xtilde=y,y=obs_[t],T=F1,Z=Z_,P=P,Q=Q,H=Hm,C=C,bUnivariate=bUnivariate,ind_non_missing=ind_non_missing_[t]) elif model_.FILTER.value == FilterAlgorithm.Durbin_Koopman.value: yt,y,residual,P,K,S,Sinv,loglk = \ dk_filter(x=yt,xtilde=y,y=obs_[t],v=residual,T=F1,Z=Z_,P=P,H=Hm,Q=Q,K=K,C=C,ind_non_missing=ind_non_missing_[t],t=t) log_likelihood += loglk return log_likelihood
[docs] def logp(p,*args): """Custom function defining log of posterioir probability.""" global count,params,param_index,model_ count += 1 #print(p) # Expected value of outcome posterior_log_likelihood = 0 try: # Get prior likelihood of model parameters prior_log_likelihood = prior_logp(p) if abs(prior_log_likelihood) > LARGE_NUMBER: return prior_log_likelihood # Update parameters parameters = np.copy(params) parameters[param_index] = p[:len(param_index)] stds = p[len(param_index):] # Get poserior likelihood of model fit to data log_likelihood = likelihood_logp(parameters,stds) # update the progress bar if PARALLEL and count%50 == 0: sys.stdout.write("\b") sys.stdout.write("=>") sys.stdout.flush() #print('parameters: ',p) #print('prior: ',prior_log_likelihood,'log_likelihood: ',log_likelihood) posterior_log_likelihood = prior_log_likelihood + log_likelihood return posterior_log_likelihood except: return -np.inf
[docs] def func_likelihood(p,*args): """Custom function defining log of likelihood.""" global count,params,param_index count += 1 # Expected value of outcome log_likelihood = 0 try: # Update parameters parameters = np.copy(params) parameters[param_index] = p[:len(param_index)] stds = p[len(param_index):] # Get poserior likelihood of model fit to data log_likelihood = likelihood_logp(parameters,stds) # update the progress bar if PARALLEL and count%50 == 0: sys.stdout.write("\b") sys.stdout.write("=>") sys.stdout.flush() return log_likelihood except: return -LARGE_NUMBER-count
[docs] def func_prior(p,*args): """Custom function defining log of prior probability.""" global count,params,param_index count += 1 #print(p) # Expected value of outcome log_likelihood = 0 try: # Get prior likelihood of model parameters log_likelihood = prior_logp(p) # update the progress bar if PARALLEL and count%50 == 0: sys.stdout.write("\b") sys.stdout.write("=>") sys.stdout.flush() return log_likelihood except: return -LARGE_NUMBER-count
[docs] def gelman_rubin(chain): """The Gelman-Rubin Test.""" ssq = np.var(chain,axis=1,ddof=1) W = np.mean(ssq,axis=0) θb = np.mean(chain,axis=1) θbb = np.mean(θb,axis=0) m = chain.shape[0] n = chain.shape[1] B = n/(m-1)*np.sum((θbb-θb)**2,axis=0) var_θ = (n-1)/n*W+1/n*B R = np.sqrt(var_θ/W) return R
[docs] def sample(model,n,obs,Qm,Hm,y0,method="emcee",parameters=None,steady_state=None, burn=10,Ndraws=310,Niter=200,ind_non_missing=None,Parallel=False, resetParameters=True,debug=True,save=False): """ Draw model parameters samples by using MCMC sampler. The variants of MC2 sampler are: * Affine Invariant Markov Chain Monte Carlo (MCMC) Ensemble sampler (emcee package). * Emcee is a Python ensemble sampling toolkit for affine-invariant MCMC. * For references on emcee package please see http://dfm.io/emcee/current/ * The algorithm is based on 2010 paper of Jonathan Goodman and Jonathan Weare’s paper, "Ensemble Samplers with Affine Invariance", https://projecteuclid.org/download/pdf_1/euclid.camcos/1513731992 * Markov Chain Monte Carlo (MCMC) sampler includes different Metropolis based sampling techniques: * Metropolis-Hastings (MH): Primary sampling method. * Adaptive-Metropolis (AM): Adapts covariance matrix at specified intervals. * Delayed-Rejection (DR): Delays rejection by sampling from a narrower distribution. Capable of n-stage delayed rejection. * Delayed Rejection Adaptive Metropolis (DRAM): DR + AM Please see https://pymcmcstat.readthedocs.io/_/downloads/en/latest/pdf/ * Markov Chain Monte Carlo (MCMC) sampler with "particles" package Please see https://pypi.org/project/particles , https://github.com/nchopin/particles , Nicolas, Chopin and Omiros, Papaspiliopoulos, 2020, "An Introduction to Sequential Monte Carlo", Springer Series in Statistics. Parameters: :param model: Model object. :type model: Model. :param n: Number of endogenous variables. :type n: int. :param obs: Measurement data. :type obs: numpy.array. :param Qm: Covariance matrix of errors of endogenous variables. :type Qm: numpy.array. :param Hm: Covariance matrix of errors of measurement variables. :type Hm: numpy.array. :param y0: Starting values of endogenous variables. :type y0: list. :param method: Sampler algorithm. :type method: str. :param parameters: List of model parameters. :type parameters: list. :param steady_state: List of steady states. :type steady_state: list. :param burn: Number of samples to discard. :type burn: int. :param Ndraws: Number of draws. :type Ndraws: int. :param Niter: Number of iterations. :type Niter: int. :param ind_non_missing: indices of non-missing observations. :type ind_non_missing: list. :param Parallel: If True runs parallel parameters sampling. :type Parallel: bool. :param resetParameters: If True resets parameters to the samples mean values. :type resetParameters: bool :param debug: If True print chain statistics information for pymcmcstat method. :type debug: bool. """ import warnings warnings.filterwarnings("ignore") global count,n_,nt,params,param_names,param_index,est_shocks_names global Z_,y0_,Qm_,obs_,data_,Hm_,ind_non_missing_,n_shocks_,shocks_ global model_,priors,pm_model,pm_names,pm_params,meas_shocks_,Pstar,PARALLEL t0 = time() n_ = n PARALLEL = Parallel model_ = model y0_ = y0; obs_= obs Qm_ = Qm; Hm_ = Hm ind_non_missing_ = ind_non_missing samples = None variables = model.calibration['variables'] n, = variables.shape shocks = meas_shocks_ = model.symbols['shocks'] n_shocks = len(shocks) n_shocks_ = n_shocks param_names = model.symbols['parameters'] params = model.calibration['parameters'] #m_params = dict(zip(param_names,params)) priors = model_.priors if 'measurement_shocks' in model.symbols: meas_shocks = meas_shocks_ = model.symbols['measurement_shocks'] else: meas_shocks = meas_shocks_ = [] # Find standard deviations of calibrated shocks est_shocks_std = [] cal = {k:v for k,v in model.priors.items() if k.startswith("std_")} for i,v in enumerate(shocks+meas_shocks): std_v = "std_"+v if std_v in cal: est_shocks_names.append(std_v) p = cal[std_v]['parameters'] est_shocks_std.append(float(p[0])) param_index = []; names = [] for i,k in enumerate(param_names): if k in priors: param_index.append(i) names.append(k) for i in range(len(est_shocks_names)): name = est_shocks_names[i] names.append(name) # Get reference to measurement function nt = len(obs) s = model.symbols['variables'] meas_variables = model.symbols['measurement_variables'] meas_vars = [x.lower() for x in meas_variables] ind = [] for mvr,x in zip(meas_variables,meas_vars): if "_meas" in x: v = mvr[:-5] elif "obs_" in x: v = mvr[4:] ind.append(s.index(v)) if 'measurement_shocks' in model.symbols: n_meas_shocks = len(model.symbols['measurement_shocks']) else: n_meas_shocks = 0 # Get reference to measurement function f_measurement = model.functions['f_measurement'] bHasAttr = hasattr(f_measurement,"py_func") mv = model.calibration['measurement_variables'] nm = len(mv) # Find measurement equations jacobian meas_params = model.calibration['measurement_parameters'] meas_var = np.zeros(n+nm+n_shocks+n_meas_shocks) if bHasAttr: meas_const,meas_jacob = f_measurement.py_func(meas_var,meas_params,order=1) else: meas_const,meas_jacob = f_measurement(meas_var,meas_params,order=1) obs -= meas_const Z = -meas_jacob[:,:n] Z_ = Z ### Get model matrices model.solved=False ls.solve(model=model,p=params,steady_state=steady_state,suppress_warnings=True) # State transition matrix F = np.copy(model.linear_model["A"]) F1 = F[:n,:n] # Array of constants C = np.copy(model.linear_model["C"][:n]) # Matrix of coefficients of shocks R = model.linear_model["R"][:n] Q = 0 for i in range(1+model.max_lead_shock-model.min_lag_shock): R1 = R[:n,i*n_shocks:(1+i)*n_shocks] Q += np.real(R1 @ Qm @ R1.T) # Set covariance matrix Nd = model.max_lead_shock - model.min_lag_shock # Get starting values of matrix P if model.PRIOR is None: Pstar = np.copy(Q) elif model.PRIOR.value == PriorAssumption.StartingValues.value: Pstar = np.copy(Q) elif model.PRIOR.value == PriorAssumption.Diffuse.value: from snowdrop.src.numeric.filters.utils import compute_Pinf_Pstar mf = np.array(np.nonzero(Z)) Pinf,Pstar = compute_Pinf_Pstar(mf=mf,T=F1,R=R,Q=Qm,N=Nd,n_shocks=n_shocks) elif model.PRIOR.value == PriorAssumption.Equilibrium.value: from snowdrop.src.numeric.solver.solvers import lyapunov_solver from snowdrop.src.numeric.solver.util import getStableUnstableRowsColumns rowStable,colStable,rowUnstable,colUnstable = getStableUnstableRowsColumns(model,T=F1,K=C) Pstar = lyapunov_solver(T=F1[np.ix_(rowStable,colStable)],R=R[rowStable],Q=Qm,N=Nd,n_shocks=n_shocks,options="discrete") elif model.PRIOR.value == PriorAssumption.Asymptotic.value: # Riccati equation from scipy.linalg import solve_discrete_are Pstar = solve_discrete_are(a=F1.T,b=Z.T,q=Q,r=Hm) ###------------------------------------------------ Sampling around optimized solution print("\nRunning Markov Chain Monte Carlo Sampling...") print(model.SAMPLING_ALGORITHM) if model.SAMPLING_ALGORITHM.value == SamplingAlgorithm.Emcee.value: # Bayesian parameter estimation with "emcee" package import emcee p = np.array(list(params[param_index]) + est_shocks_std) ndim = len(p) nwalkers = 2*ndim pos = p * (1+0.5*np.random.rand(nwalkers,ndim)) if PARALLEL: n_threads = num_cores else: n_threads = 1 # Set sampler if emcee.__version__ < "3": sampler = emcee.EnsembleSampler(nwalkers=nwalkers,dim=ndim,lnpostfn=logp,threads=n_threads) else: sampler = emcee.EnsembleSampler(nwalkers=nwalkers,ndim=ndim,log_prob_fn=logp,threads=n_threads) # Run the MCMC for N steps sampler.run_mcmc(pos, Ndraws, progress=True) chain = sampler.chain # Discard the initial "burn" steps samples = chain[:, burn:, :] samples = samples.reshape((-1, ndim)) panel = np.mean(chain,axis=0) panel2 = chain.reshape((-1, ndim)) params_last = chain[-1,-1,:] if debug: print("\n\nAcceptance fraction: {:.1f}%".format(100*np.mean(sampler.acceptance_fraction))) print("\n\tμ\t\t\tσ\t\t\tGelman-Rubin") print("========================================") print("\t{:.2f}\t\t{:.2f}\t\t{:.2f}".format( chain.reshape(-1, chain.shape[-1]).mean(axis=0)[0], chain.reshape(-1, chain.shape[-1]).std(axis=0)[0], gelman_rubin(chain)[0])) from snowdrop.src.graphs.util import plot_chain,plot_pairwise_correlation plot_chain(path_to_dir=path_to_dir,chains=panel,names=names,title="Convergence",save=save) plot_pairwise_correlation(path_to_dir=path_to_dir,chains=panel2,names=names,title="Pairwise Correlation",save=save) elif model.SAMPLING_ALGORITHM.value in [SamplingAlgorithm.Particle_pmmh.value, SamplingAlgorithm.Particle_smc.value, SamplingAlgorithm.Particle_gibbs.value]: # Bayesian parameter estimation with "particles" package from snowdrop.src.utils.distributions import StructDist #import particles.distributions as dists # Distributions from particles import mcmc # MCMC algorithms (PMMH, Particle Gibbs, etc.) from particles import smc_samplers as ssp #from particles import collectors as col from collections import OrderedDict from snowdrop.src.utils.distributions import getParticleParameter from snowdrop.src.numeric.filters.particle import StateSpaceModel from snowdrop.src.numeric.filters.particle import ParticleGibbs #from statsmodels.stats.correlation_tools import cov_nearest as nearestPositiveDefinite non_empty = [bool(x) for x in ind_non_missing_] data = data_ = obs[non_empty] ### Get matrices # State transition matrix. #T_ = model.linear_model["A"][:n,:n] # Array of constants. #C_ = model.linear_model["C"][:n] # Matrix of coefficients of shocks R = model.linear_model["R"][:n] Q_ = 0 for i in range(1+model.max_lead_shock-model.min_lag_shock): R1 = R[:,i*n_shocks:(1+i)*n_shocks] Q_ += np.real(R1 @ Qm @ R1.T) Q_ = 0.5*(Q_+Q_.T) # Initialize parameter's array prior_names = []; prior_dict = OrderedDict(); val = [] for i,index in enumerate(param_index): name = param_names[index] prior = model_.priors[name] pars = prior['parameters'] distr = prior['distribution'] p = getParticleParameter(distr,pars) prior_names.append(name) prior_dict[name] = p val.append(float(pars[0])) for i in range(len(est_shocks_names)): name = est_shocks_names[i] prior = model_.priors[name] pars = prior['parameters'] distr = prior['distribution'] p = getParticleParameter(distr,pars) prior_dict[name] = p my_prior = StructDist(prior_dict) # Initial values of parameters theta0 = np.zeros((1), dtype=my_prior.dtype) for v,name in zip(val,names): theta0[name] = v for i in range(len(est_shocks_names)): theta0[est_shocks_names[i]] = est_shocks_std[i] # Initial number of particles Nparticles = 100 # Initialize state-space model # ss_model = StateSpaceModel(x=y0,y=obs,theta0=theta0,data=data,T=T_,R=R,Z=Z,P=Q_,Q=Q_,Qm=Qm,H=Hm,C=C_,n_shocks=n_shocks,verbose=True) # Simulate from the model 'nt' data points # states, y_data = ss_model.simulate(nt) if model.SAMPLING_ALGORITHM.value == SamplingAlgorithm.Particle_gibbs.value: # Particle Gibbs sampler for a state-space model pf = ParticleGibbs(niter=Niter,Nx=Nparticles,ssm_cls=StateSpaceModel, theta0=theta0,prior=my_prior,data=data) pf.run() theta = pf.chain.theta elif model.SAMPLING_ALGORITHM.value == SamplingAlgorithm.Particle_pmmh.value: # Particle Marginal Metropolis Hastings algorithm pf = mcmc.PMMH(niter=Niter,Nx=Ndraws,ssm_cls=StateSpaceModel, prior=my_prior,data=data,theta0=theta0) pf.run() theta = pf.chain.theta elif model.SAMPLING_ALGORITHM.value == SamplingAlgorithm.Particle_smc.value: # Sequential Quasi-Monte Carlo algorithm import particles fk_smc = ssp.SMC2(ssm_cls=StateSpaceModel,prior=my_prior,data=data, init_Nx=Nparticles,wastefree=True) collect = None #[col.Moments()] pf = particles.SMC(fk=fk_smc,qmc=False,N=Nparticles,store_history=False, collect=collect,verbose=False) pf.run() # for i in range(len(data)-1): # next(pf) theta = pf.X.theta samples = [] for name in names: z = theta[name] if len(z) > burn: y = z[burn:] # discard the first "burn" iterations else: y = z samples.append(y) samples = np.array(samples).T params_last = samples[-1] if debug: from snowdrop.src.graphs.util import plot_chain,plot_pairwise_correlation plot_chain(path_to_dir=path_to_dir,chains=samples,names=names,title="Convergence",save=save) plot_pairwise_correlation(path_to_dir=path_to_dir,chains=samples,names=names,title="Pairwise Correlation",save=save) elif model.SAMPLING_ALGORITHM.value in [SamplingAlgorithm.Pymcmcstat_mh.value, SamplingAlgorithm.Pymcmcstat_am.value, SamplingAlgorithm.Pymcmcstat_dr.value, SamplingAlgorithm.Pymcmcstat_dram.value]: # Bayesian parameter estimation with "pymcmcstat" package import tempfile from datetime import datetime import pymcmcstat as mcstat from pymcmcstat.MCMC import MCMC from pymcmcstat.ParallelMCMC import ParallelMCMC from pymcmcstat.chain import ChainProcessing as CP from pymcmcstat.chain import ChainStatistics as CS # Define a directory name based on the date/time to ensure we do not overwrite any results. datestr = datetime.now().strftime('%Y%m%d_%H%M%S') savedir = tempfile.gettempdir() + os.sep + 'resources' + os.sep + str('{}_{}'.format(datestr, 'parallel_chains')) if model.SAMPLING_ALGORITHM.value == SamplingAlgorithm.Pymcmcstat_mh.value: sampling_method = "mh" elif model.SAMPLING_ALGORITHM.value == SamplingAlgorithm.Pymcmcstat_am.value: sampling_method = "am" elif model.SAMPLING_ALGORITHM.value == SamplingAlgorithm.Pymcmcstat_dr.value: sampling_method = "dr" elif model.SAMPLING_ALGORITHM.value == SamplingAlgorithm.Pymcmcstat_dram.value: sampling_method = "dram" # Initialize MCMC object mc2 = MCMC() # initialize data structure mc2.data.add_data_set(x=[], y=[]) mc2.simulation_options.define_simulation_options( nsimu=Ndraws,updatesigma=True,method=sampling_method, savedir=savedir,savesize=1000,save_to_json=True, verbosity=0,waitbar=not PARALLEL,save_lightly=True,save_to_bin=True) mc2.model_settings.define_model_settings(sos_function=func_likelihood,prior_function=func_prior,sigma2=[1]) # Initialize parameter's array iv = [] for i,index in enumerate(param_index): name = param_names[index] prior = model_.priors[name] pars = prior['parameters'] lower_bound = float(pars[1]) upper_bound = float(pars[2]) x = pars[0] #x = m_params[name] iv.append(x) #print(name,x,lower_bound,upper_bound) mc2.parameters.add_model_parameter(name=name,theta0=x,minimum=lower_bound,maximum=upper_bound,sample=1) # Initialize standard deviation's array for name in est_shocks_names: prior = model_.priors[name] pars = prior['parameters'] lower_bound = float(pars[1]) upper_bound = float(pars[2]) x = pars[0] iv.append(x) #print(name,x,lower_bound,upper_bound) mc2.parameters.add_model_parameter(name=name,theta0=x,minimum=lower_bound,maximum=upper_bound,sample=1) if PARALLEL: # # Setup Parallel Simulation and Define Initial Values # - You can specify the number of chains to be generated (`num_chain`) # and the number of cores to use (`num_cores`). # setup parallel MCMC parMC = ParallelMCMC() parMC.setup_parallel_simulation(mcset=mc2,initial_values=np.array([iv]*num_chains), num_chain=num_chains,num_cores=num_cores) try: # Run Parallel Simulations. parMC.run_parallel_simulation() # Display Chain Statistics. #parMC.display_individual_chain_statistics() # Processing Chains # To load the results we must specify the name of the top level directory # where the parallel results are stored. pres = CP.load_parallel_simulation_results(savedir) combined_chains,index = CP.generate_combined_chain_with_index(pres) chains = CP.generate_chain_list(pres,burnin_percentage=int(100.*burn/Ndraws)) # Aggregate chains samples = np.concatenate(chains) params_last = samples[-1] if debug: # Gelman-Rubin Diagnostic # We generated the list of chains which are required for the Gelman-Rubin diagnostic. # The output will be a dictionary where the keys are the parameter names. # Each parameter references another dictionary which contains print(" - `R`, Potential Scale Reduction Factor (PSRF)") print(" - `B`, Between Sequence Variance") print(" - `W`, Within Sequence Variance") print(" - `V`, Mixture-of-Sequences Variance") print(" - `neff`, Effective Number of Samples") CS.gelman_rubin(chains=chains, display=True) from snowdrop.src.graphs.util import plot_chain,plot_pairwise_correlation plot_chain(path_to_dir=path_to_dir,chains=chains[-1],names=names,title="Convergence",save=save) plot_pairwise_correlation(path_to_dir=path_to_dir,chains=chains[-1],names=names,title="Pairwise Correlation",save=save) except: samples = [] finally: # Remove temporary directory if os.path.exists(savedir): import shutil shutil.rmtree(savedir) else: try: # Run Sequential Simulations. mc2.run_simulation() # Extract results. results = mc2.simulation_results.results chain = results['chain'] samples = chain[burn:] params_last = chain[-1] if debug: # Display chain statistics mc2.chainstats(samples,results) # Generate mcmc plots mcpl = mcstat.mcmcplot # initialize plotting methods mcpl.plot_chain_panel(samples,names) mcpl.plot_density_panel(samples,names) mcpl.plot_pairwise_correlation_panel(samples,names) except Exception as ex: print(ex) finally: # Remove temporary directory if os.path.exists(savedir): import shutil shutil.rmtree(savedir) else: raise Exception(f"Undefined sampling algorithm: {model.SAMPLING_ALGORITHM.name}") elapsed = time() - t0 elapsed = np.round(elapsed,2) # Mean, median and standard deviation of sample params_mean = np.mean(samples, axis=0) params_median = np.median(samples, axis=0) params_std = np.std(samples, axis=0) header = ['Name','Start Val.','Median','Std.','Last','Shape','L.B.','U.B.'] pTable = PrettyTable(header) pTable.float_format = '5.3' # Estimated parameters p = params[param_index] for i,name in enumerate(names): prior = model.priors[name] distr = prior['distribution'] pars = prior['parameters'] par = p[i] if i < len(p) else pars[0] row = [name,par,params_median[i],params_std[i],params_last[i],distr.replace("_pdf","").replace("_hp",""),pars[1],pars[2]] pTable.add_row(row) cprint("\nEnsemble Estimation","green") print(pTable) # Save parameters all_parameters = np.copy(parameters) if resetParameters: all_parameters[param_index] = params_mean #all_parameters[param_index] = params_last sys.stdout.write("\n\n") # this ends the progress bar model.solved = False print("\nSampling elapsed time: {0} (seconds)".format(elapsed)) return all_parameters,names,samples