diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index ec120c5e73fa1491f8984b5ed8d2f02bc3a79eb7..8c52958d320c8929efc684ddf26f0769597ddd49 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -41,7 +41,7 @@ class BayesInference:
self.PCEStd = {}
self.Discrepancy = None
self.ReqOutputType = None
- self.NrofSamples = None
+ self.NrofSamples = 100000
self.Samples = None
self.Samplesu = None
self.MCMCnSteps = None
@@ -104,9 +104,9 @@ class BayesInference:
meanPCEOutputs = {}
stdPCEOutputs = {}
univ_p_val = PCEModel.univ_basis_vals(Samples)
+
for Outkey, ValuesDict in PCEModel.CoeffsDict.items():
-
PCEOutputs_mean = np.zeros((len(Samples), len(ValuesDict)))
PCEOutputs_std = np.zeros((len(Samples), len(ValuesDict)))
idx = 0
@@ -115,14 +115,15 @@ class BayesInference:
PSI_Val = PCEModel.PCE_create_Psi(PolynomialDegrees, univ_p_val)
# Perdiction
- try: #with error bar
+ if PCEModel.RegMethod.lower() == 'ols': #without error bar
+ coeffs = PCEModel.CoeffsDict[Outkey][Inkey]
+ PCEOutputs_mean[:, idx] = np.dot(PSI_Val, coeffs)
+
+ else: #with error bar
clf_poly = PCEModel.clf_poly[Outkey][Inkey]
y_mean, y_std = clf_poly.predict(PSI_Val, return_std=True)
PCEOutputs_mean[:, idx] = y_mean
PCEOutputs_std[:, idx] = y_std
- except: #without error bar
- coeffs = PCEModel.CoeffsDict[Outkey][Inkey]
- PCEOutputs_mean[:, idx] = np.dot(PSI_Val, coeffs)
if PCEModel.metaModel == 'PCEKriging':
gp = PCEModel.gp_poly[Outkey][Inkey]
@@ -132,13 +133,11 @@ class BayesInference:
idx += 1
-
if PCEModel.index is None:
if PCEModel.DimRedMethod.lower() == 'pca':
PCA = PCEModel.pca[Outkey]
meanPCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_)
stdPCEOutputs[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))
- #np.dot(PCEOutputs_std, PCA.components_)
else:
meanPCEOutputs[Outkey] = PCEOutputs_mean
stdPCEOutputs[Outkey] = PCEOutputs_std
@@ -173,8 +172,7 @@ class BayesInference:
Model = self.PCEModel.ModelObj
if Samples is None:
- PCEModel = self.PCEModel
- Samples = self.get_Sample(PCEModel)
+ Samples = self.get_Sample()
self.Samples = Samples
else:
Samples = Samples
@@ -547,10 +545,6 @@ class BayesInference:
self.Samples = SamplesAllParameters
self.NrofSamples = self.Samples.shape[0]
- NrofSamples = self.NrofSamples
- else:
- NrofSamples = self.NrofSamples
-
# ---------- Preparation of observation data ----------
# Read observation data and perturb it if requested
@@ -619,57 +613,57 @@ class BayesInference:
Data[itrIdx] = data.flatten('F')
- # -----------------------------------------------------
- # ----- Loop over the perturbed observation data ------
- # -----------------------------------------------------
- # Initilize arrays
- Likelihoods = np.zeros((NrofSamples,BootstrapItrNr))
- BME_Corr = np.zeros((BootstrapItrNr))
- BME = np.zeros((BootstrapItrNr))
-
- for itrIdx, data in tqdm(enumerate(Data), ascii=True, desc ="Boostraping the BME calculations"):
+ # -----------------------------------------------------
+ # ----- Loop over the perturbed observation data ------
+ # -----------------------------------------------------
+ # Initilize arrays
+ Likelihoods = np.zeros((self.NrofSamples,BootstrapItrNr))
+ BME_Corr = np.zeros((BootstrapItrNr))
+ BME = np.zeros((BootstrapItrNr))
- # ---------------- Likelihood calculation ----------------
- if self.emulator:
- # Evaluate the PCEModel
- self.meanPCEPriorPred, self.stdPCEPriorPred = self.eval_PCEmodel(Samples=self.Samples)
+ for itrIdx, data in tqdm(enumerate(Data), ascii=True, desc ="Boostraping the BME calculations"):
- # unknown sigma2
- if optSigma == 'C':
- Likelihoods[:,itrIdx] = self.normpdfSigma2(self.meanPCEPriorPred, data, TotalSigma2)
- else:
- # known sigma2
- Likelihoods[:,itrIdx] = self.normpdf(self.meanPCEPriorPred, data, TotalSigma2)
-
- else:
- self.ModelPriorPred = self.eval_Model(Samples=self.Samples)
- # unknown sigma2
- if optSigma == 'C':
- Likelihoods[:,itrIdx] = self.normpdfSigma2(self.ModelPriorPred, data, TotalSigma2)
+ # ---------------- Likelihood calculation ----------------
+ if self.emulator:
+ # Evaluate the PCEModel
+ self.meanPCEPriorPred, self.stdPCEPriorPred = self.eval_PCEmodel(Samples=self.Samples)
+
+ # unknown sigma2
+ if optSigma == 'C':
+ Likelihoods[:,itrIdx] = self.normpdfSigma2(self.meanPCEPriorPred, data, TotalSigma2)
+ else:
+ # known sigma2
+ Likelihoods[:,itrIdx] = self.normpdf(self.meanPCEPriorPred, data, TotalSigma2)
+
else:
- # known sigma2
- Likelihoods[:,itrIdx] = self.normpdf(self.ModelPriorPred, data, TotalSigma2)
-
- # ---------------- BME Calculations ----------------
- # BME (log)
- BME[itrIdx] = np.log(np.nanmean(np.exp(Likelihoods[:,itrIdx])))
+ self.ModelPriorPred = self.eval_Model(Samples=self.Samples)
+ # unknown sigma2
+ if optSigma == 'C':
+ Likelihoods[:,itrIdx] = self.normpdfSigma2(self.ModelPriorPred, data, TotalSigma2)
+ else:
+ # known sigma2
+ Likelihoods[:,itrIdx] = self.normpdf(self.ModelPriorPred, data, TotalSigma2)
+
+ # ---------------- BME Calculations ----------------
+ # BME (log)
+ BME[itrIdx] = np.log(np.nanmean(np.exp(Likelihoods[:,itrIdx])))
+
+ # BME correction when using Emulator
+ # if self.emulator:
+ # BME_Corr[itrIdx] = self.BME_Corr_Weight(data, sigma2)
+
+ # ---------------- Store BME, Likelihoods for all ----------------
+ # Likelihoods (Size: NrofSamples,BootstrapItrNr)
+ self.Likelihoods = Likelihoods
- # BME correction when using Emulator
- # if self.emulator:
- # BME_Corr[itrIdx] = self.BME_Corr_Weight(data, sigma2)
-
- # ---------------- Store BME, Likelihoods for all ----------------
- # Likelihoods (Size: NrofSamples,BootstrapItrNr)
- self.Likelihoods = Likelihoods
-
- # BME (log) (Size: 1,BootstrapItrNr)
- self.BME = BME
-
- # TODO: BMECorrFactor (log) (Size: 1,BootstrapItrNr)
- #if self.emulator: self.BMECorrFactor = BME_Corr
-
- # BME = BME + BMECorrFactor
- if self.emulator: self.BME = self.BME #+ self.BMECorrFactor
+ # BME (log) (Size: 1,BootstrapItrNr)
+ self.BME = BME
+
+ # TODO: BMECorrFactor (log) (Size: 1,BootstrapItrNr)
+ #if self.emulator: self.BMECorrFactor = BME_Corr
+
+ # BME = BME + BMECorrFactor
+ if self.emulator: self.BME = self.BME #+ self.BMECorrFactor
# ---------------- Parameter Bayesian inference ----------------
if self.SamplingMethod.lower() == 'mcmc':
@@ -686,7 +680,7 @@ class BayesInference:
# -------- Find MAP and run PCEModel and origModel --------
# Compute the MAP
- if self.MAP =='Mean':
+ if self.MAP.lower() =='mean':
MAP_theta = self.Posterior_df.to_numpy().mean(axis=0).reshape((1,NofPa))
else:
MAP_theta = stats.mode(self.Posterior_df.to_numpy(),axis=0)[0]
@@ -715,24 +709,25 @@ class BayesInference:
# -------- Posteior parameters --------
parNames = [self.PCEModel.Inputs.Marginals[i].Name for i in range(NofPa)]
figPosterior = corner.corner(self.Posterior_df.to_numpy(), labels=parNames,
- show_titles=True, title_fmt=self.Corner_title_fmt,
- title_kwargs={"fontsize": 14})
-
- # Extract the axes
- axes = np.array(figPosterior.axes).reshape((NofPa, NofPa))
-
- # Loop over the diagonal
- for i in range(NofPa):
- ax = axes[i, i]
- ax.axvline(MAP_theta[0,i], ls='--', color="r")
-
- # Loop over the histograms
- for yi in range(NofPa):
- for xi in range(yi):
- ax = axes[yi, xi]
- ax.axvline(MAP_theta[0,xi], ls='--', color="r")
- ax.axhline(MAP_theta[0,yi], ls='--', color="r")
- ax.plot(MAP_theta[0,xi], MAP_theta[0,yi], "sr")
+ quantiles=[0.15, 0.5, 0.85],show_titles=True,
+ title_fmt=self.Corner_title_fmt,
+ title_kwargs={"fontsize": 12})
+
+ # # Extract the axes
+ # axes = np.array(figPosterior.axes).reshape((NofPa, NofPa))
+
+ # # Loop over the diagonal
+ # for i in range(NofPa):
+ # ax = axes[i, i]
+ # ax.axvline(MAP_theta[0,i], ls='--', color="r")
+
+ # # Loop over the histograms
+ # for yi in range(NofPa):
+ # for xi in range(yi):
+ # ax = axes[yi, xi]
+ # ax.axvline(MAP_theta[0,xi], ls='--', color="r")
+ # ax.axhline(MAP_theta[0,yi], ls='--', color="r")
+ # ax.plot(MAP_theta[0,xi], MAP_theta[0,yi], "sr")
figPosterior.set_size_inches((24,16))
diff --git a/BayesValidRox/BayesInference/MCMC.py b/BayesValidRox/BayesInference/MCMC.py
index 74a6190c87272534ec305e1c8e245e7823fafc7c..42e875bca9ab4ed44d0ecd60562911dae44e737e 100755
--- a/BayesValidRox/BayesInference/MCMC.py
+++ b/BayesValidRox/BayesInference/MCMC.py
@@ -13,7 +13,6 @@ from multiprocessing import Pool
import scipy.stats as st
from scipy.linalg import cholesky as chol
import warnings
-import time
class MCMC():
"""
@@ -49,7 +48,11 @@ class MCMC():
# set theta near the maximum likelihood, with
np.random.seed(0)
if self.initsamples is None:
- initsamples = priorDist.sample(self.nwalkers).T
+ try:
+ initsamples = priorDist.sample(self.nwalkers).T
+ except:
+ # when aPCE selected - gaussian kernel distribution
+ initsamples = priorDist.resample(self.nwalkers).T
else:
# Pick samples based on a uniform dist between main and max of each dim
initsamples = np.zeros((self.nwalkers, ndim))
@@ -109,8 +112,8 @@ class MCMC():
for i in range(nsamples):
if not self.check_ranges(theta[i], params_range):
- logprior[i] = priorDist.pdf(theta[i])
-
+ logprior[i] = np.log(priorDist.pdf(theta[i]))
+
if nsamples == 1:
return logprior[0]
else:
@@ -129,10 +132,9 @@ class MCMC():
else:
# Evaluate the origModel
meanPCEPriorPred = BayesOpts.eval_Model(theta)
-
+
return BayesOpts.normpdf(meanPCEPriorPred, Observation, TotalSigma2)
-
-
+
def log_posterior(self, theta, Observation, TotalSigma2):
"""
Computes the posterior likelihood for the given parameterset.
@@ -152,6 +154,7 @@ class MCMC():
DESCRIPTION.
"""
+
ndimTheta = theta.ndim
nsamples = 1 if ndimTheta == 1 else theta.shape[0]
diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index e9eae4e6c5fa37239b03de814f7a6691597ba4d8..776f7520baea9aa6b2ec427922c60bdfa300e89f 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -858,13 +858,12 @@ class PostProcessing:
# Initialize the cell to store the names of the variables of the Sobol Idx.
TotalVariance = np.zeros((NofMeasurements))
-
# Loop over all measurement points and calculate sobol indices
for PointIdx in range(NofMeasurements):
# Extract the basis indices (alpha) and coefficients
PolynomialDegree = PCEModel.BasisDict[Output]['y_%s'%(PointIdx+1)]
-
+
try:
clf_poly = PCEModel.clf_poly[Output]['y_%s'%(PointIdx+1)]
PCECoeffs = clf_poly.coef_
@@ -873,7 +872,7 @@ class PostProcessing:
# Compute total variance
TotalVariance[PointIdx] = np.sum(np.square(PCECoeffs[1:]))
-
+
nzidx = np.where(PCECoeffs != 0)[0]
# Set all the Sobol indices equal to zero if the presence of a null output.
if len(nzidx) == 0:
@@ -902,7 +901,7 @@ class PostProcessing:
sum_ind = nzidx[idx]
total_sobol_array[ParIdx,PointIdx] = np.sum(np.square(PCECoeffs[sum_ind])) / TotalVariance[PointIdx]
-
+
# ----- if PCA selected: Compute covariance # -----
if PCEModel.DimRedMethod.lower() == 'pca':
cov_Z_p_q = np.zeros((NofPa))
@@ -926,8 +925,7 @@ class PostProcessing:
cov_Z_p_q[ParIdx] = np.sum(PCECoeffs[sum_ind] * nextPCECoeffs[sum_ind])
except:
cov_Z_p_q[ParIdx] = 0.0
-
-
+
# Compute the sobol indices according to Ref. 2
if PCEModel.DimRedMethod.lower() == 'pca':
ntMeasurements = PCEModel.ExpDesign.Y[Output].shape[1]
diff --git a/BayesValidRox/surrogate_models/ExpDesigns.py b/BayesValidRox/surrogate_models/ExpDesigns.py
index 53a21c9c7fd4473cdfd1aacf409adc5527d771b3..3707a5485f9a510bbd0379e4790daaec46f9912a 100644
--- a/BayesValidRox/surrogate_models/ExpDesigns.py
+++ b/BayesValidRox/surrogate_models/ExpDesigns.py
@@ -10,6 +10,7 @@ import numpy as np
import math
import itertools
import chaospy
+import scipy.stats as st
from tqdm import tqdm
from .aPoly_Construction import aPoly_Construction
@@ -64,7 +65,7 @@ class ExpDesigns:
NofPa = len(Inputs.Marginals)
self.SamplingMethod = SamplingMethod
NrSamples = int(NrSamples)
-
+ self.arbitrary = False if len(Inputs.Marginals[0].InputValues) == 0 else True
# Get the bounds if InputValues are directly defined by user:
if Inputs.Marginals[0].Parameters is None:
for i in range(NofPa):
@@ -75,7 +76,7 @@ class ExpDesigns:
self.JDist, self.Polytypes = self.DistConstructor()
## Sample the distribution of parameters
- if len(Inputs.Marginals[0].InputValues) == 0 and self.SamplingMethod != 'user':
+ if not self.arbitrary and self.SamplingMethod != 'user':
## Case I = polytype not arbitrary
# Execute initialization to get the boundtuples
@@ -88,7 +89,7 @@ class ExpDesigns:
else:
# Case II:
# polytype arbitrary or Input values are directly prescribed by the user.
-
+
# Generate the samples based on requested method
if self.SamplingMethod == 'user':
X = self.X
@@ -105,7 +106,10 @@ class ExpDesigns:
else:
raise Exception('Unfortunately, for the prescribed inputs, the samepling method '
'%s is not valid.'%self.SamplingMethod)
-
+
+ # Naive approach: Fit a gaussian kernel to the provided data
+ self.JDist = st.gaussian_kde(self.Input_distributions)
+
return X
def DistConstructor(self):