diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index 1b291e03f8983316af2a2c06f78d08eae01f59b5..70404e348784e2771d7dcd12a5a7a5ae670dacf6 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -10,6 +10,7 @@ import os, sys
import pandas as pd
from tqdm import tqdm
from scipy import stats
+import scipy.linalg as spla
import seaborn as sns
import emcee
import corner
@@ -51,6 +52,7 @@ class BayesInference:
self.MCMCnwalkers = None
self.MCMCinitSamples = None
self.MCMCmoves = None
+ self.MCMCverbose = False
self.ModelOutputs = []
self.meanPCEPriorPred = []
self.stdPCEPriorPred = []
@@ -172,7 +174,7 @@ class BayesInference:
covMatrix += covBias
# Add the std of the PCE is chosen as emulator.
- if self.emulator and self.PCEModel.metaModel.lower() != 'pcekriging':
+ if self.emulator:
covMatrix_PCE = np.zeros((NofMeasurements, NofMeasurements), float)
stdPCE = np.empty((SampleSize,0))
for outputType in OutputType:
@@ -191,13 +193,19 @@ class BayesInference:
covMatrix = covMatrix[self.selectedIndices,self.selectedIndices]
# Compute loglikelihood
- try:
- logL = multivariate_normal.logpdf(TotalOutputs, mean=Data, cov=covMatrix)
- except:
- logL = -np.inf
-
- return logL
-
+ if TotalOutputs.shape[0] == 1:
+ return self._logpdf(TotalOutputs[0], Data, covMatrix)
+ else:
+ return np.array([self._logpdf(output, Data, covMatrix) for output in TotalOutputs])
+
+ def _logpdf(self, x, mean, cov):
+ n = len(mean)
+ L = spla.cholesky(cov, lower=True)
+ beta = np.sum(np.log(np.diag(L)))
+ dev = x - mean
+ alpha = dev.dot(spla.cho_solve((L, True), dev))
+ return -0.5 * alpha - beta - n / 2. * np.log(2 * np.pi)
+
#--------------------------------------------------------------------------------------------------------
def normpdfSigma2(self, Outputs, Data, Sigma2):
@@ -655,7 +663,7 @@ class BayesInference:
nsteps = 100000 if self.MCMCnSteps is None else self.MCMCnSteps
nwalkers = 50*self.PCEModel.NofPa if self.MCMCnwalkers is None else self.MCMCnwalkers
multiprocessing = False if self.MultiProcessMCMC is None else self.MultiProcessMCMC
- MCMC_ = MCMC(self, initsamples=initsamples, nwalkers=nwalkers,
+ MCMC_ = MCMC(self, initsamples=initsamples, nwalkers=nwalkers, verbose=self.MCMCverbose,
nsteps = nsteps, moves=self.MCMCmoves, multiprocessing=multiprocessing)
self.Posterior_df = MCMC_.run_sampler(Data[0], TotalSigma2)
diff --git a/BayesValidRox/BayesInference/MCMC.py b/BayesValidRox/BayesInference/MCMC.py
index 96f4afee65bd0653bfac4ec0c1fe869bb48fd4ae..7dc0d88044b83499734f8e9b91c36330d8b49444 100755
--- a/BayesValidRox/BayesInference/MCMC.py
+++ b/BayesValidRox/BayesInference/MCMC.py
@@ -27,7 +27,7 @@ class MCMC():
ntemps = 20
"""
def __init__(self, BayesOpts, initsamples = None, nwalkers = 100, nburn = 100,
- nsteps = 100000, moves = None, multiprocessing=False, verbose = False):
+ nsteps = 100000, moves = emcee.moves.KDEMove(), multiprocessing=False, verbose = False):
self.BayesOpts = BayesOpts
self.initsamples = initsamples
@@ -159,7 +159,7 @@ class MCMC():
# always print the current mean acceptance fraction
print("\nStep: {}".format(sampler.iteration))
- print("Mean acceptance fraction: {0:.3f}".format(np.mean(sampler.acceptance_fraction)))
+ # print("Mean acceptance fraction: {0:.3f}".format(np.mean(sampler.acceptance_fraction)))
# compute the autocorrelation time so far
# using tol=0 means that we'll always get an estimate even if it isn't trustworthy
@@ -173,8 +173,6 @@ class MCMC():
# check convergence
converged = np.all(tau * 100 < sampler.iteration)
converged &= np.all(np.abs(tauold - tau) / tau < 0.01)
- # converged &= np.all(sampler.acceptance_fraction < 0.5)
- # converged &= np.all(sampler.acceptance_fraction > 0.234)
if converged:
break
tauold = tau
@@ -198,10 +196,10 @@ class MCMC():
print("thin: {0}".format(thin))
print("burn-in: {0}".format(burnin))
print("flat chain shape: {0}".format(finalsamples.shape))
- print("Mean acceptance fraction*: {0:.3f}".format(np.nanmean(sampler.acceptance_fraction)))
- print("Gelman-Rubin Test**: ", np.round(self.gelman_rubin(sampler.chain[:,burnin:]),3))
- print("\n* This value must lay between 0.234 and 0.5.")
- print("** These values must be smaller than 1.1.")
+ # print("Mean acceptance fraction*: {0:.3f}".format(np.nanmean(sampler.acceptance_fraction)))
+ print("Gelman-Rubin Test*: ", np.round(self.gelman_rubin(sampler.chain[:,burnin:]),3))
+ # print("\n* This value must lay between 0.234 and 0.5.")
+ print("* These values must be smaller than 1.1.")
print('-'*50)
print("\n>>>> Bayesian inference with MCMC for {0} successfully completed. <<<<<<\n".format(self.BayesOpts.Name))
@@ -215,42 +213,43 @@ class MCMC():
# create a PdfPages object
- pdf = PdfPages(OutputDir+'/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, '-')
+ if self.verbose:
+ pdf = PdfPages(OutputDir+'/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 ' + parNames[parIdx])
+ main_ax.set_xlabel('step number')
- # histogram on the attached axes
- y_hist.hist(samples[burnin:], 40, histtype='stepfilled',
- orientation='horizontal', color='gray')
+ # save the current figure
+ pdf.savefig(fig, bbox_inches='tight')
- main_ax.set_ylim(params_range[parIdx])
- main_ax.set_title('traceplot for ' + parNames[parIdx])
- main_ax.set_xlabel('step number')
+ # Destroy the current plot
+ plt.clf()
- # save the current figure
- pdf.savefig(fig, bbox_inches='tight')
-
- # Destroy the current plot
- plt.clf()
-
- pdf.close()
+ pdf.close()
# plot development of autocorrelation estimate
if not self.mp:
fig1 = plt.figure()
- steps = autocorreverynsteps*np.arange(1, autocorrIdx+1)
+ steps = 100*np.arange(1, autocorrIdx+1)
taus = autocorr[:autocorrIdx]
- plt.plot(steps, steps / float(autocorreverynsteps), "--k")
+ plt.plot(steps, steps / 100.0, "--k")
plt.plot(steps, taus)
plt.xlim(0, steps.max())
plt.ylim(0, np.nanmax(taus) + 0.1*(np.nanmax(taus) - np.nanmin(taus)))