diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index b6697b92af8e72921d39b37167139f1d62a707f5..ec120c5e73fa1491f8984b5ed8d2f02bc3a79eb7 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -450,10 +450,7 @@ class BayesInference:
for name in self.Discrepancy.Name:
InputNames.append(name)
- self.Posterior_df = pd.DataFrame(acceptedSamples, columns=InputNames)
-
-
- return self.Posterior_df
+ return pd.DataFrame(acceptedSamples, columns=InputNames)
#--------------------------------------------------------------------------------------------------------
def PosteriorPredictive(self):
@@ -625,59 +622,57 @@ class BayesInference:
# -----------------------------------------------------
# ----- Loop over the perturbed observation data ------
# -----------------------------------------------------
- if self.SamplingMethod == 'rejection':
- # Initilize arrays
- Likelihoods = np.zeros((NrofSamples,BootstrapItrNr))
- BME_Corr = np.zeros((BootstrapItrNr))
- BME = np.zeros((BootstrapItrNr))
+ # 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"):
- for itrIdx, data in tqdm(enumerate(Data), ascii=True, desc ="Boostraping the BME calculations"):
+ # ---------------- Likelihood calculation ----------------
+ if self.emulator:
+ # Evaluate the PCEModel
+ self.meanPCEPriorPred, self.stdPCEPriorPred = self.eval_PCEmodel(Samples=self.Samples)
- # ---------------- 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)
-
+ # unknown sigma2
+ if optSigma == 'C':
+ Likelihoods[:,itrIdx] = self.normpdfSigma2(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)
- 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 (log) (Size: 1,BootstrapItrNr)
- self.BME = BME
+ # known sigma2
+ Likelihoods[:,itrIdx] = self.normpdf(self.meanPCEPriorPred, data, TotalSigma2)
- # 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
+ else:
+ 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)
- # ---------------- Rejection Sampling ----------------
- self.Posterior_df = self.Rejection_Sampling()
- else:
+ # ---------------- 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 (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':
# TODO: MCMC
initsamples = None if self.Samples is None else self.Samples
nsteps = 1000 if self.MCMCnSteps is None else self.MCMCnSteps
@@ -685,7 +680,9 @@ class BayesInference:
MCMC_ = MCMC(self, initsamples=initsamples, nwalkers=2*self.PCEModel.NofPa,
nsteps = nsteps, multiprocessing=multiprocessing)
self.Posterior_df = MCMC_.run_sampler(Data[0], TotalSigma2)
-
+
+ else: # Rejection sampling
+ self.Posterior_df = self.Rejection_Sampling()
# -------- Find MAP and run PCEModel and origModel --------
# Compute the MAP
diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index 613fea138b228c2e193bc669e1875d7b3a92c507..e9eae4e6c5fa37239b03de814f7a6691597ba4d8 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -185,7 +185,7 @@ class PostProcessing:
if PCEModel.index is None:
if PCEModel.DimRedMethod.lower() == 'pca':
PCA = PCEModel.pca[Outkey]
- self.PCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_)
+ self.PCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_)
self.PCEOutputs_std[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))
else:
self.PCEOutputs[Outkey] = PCEOutputs_mean
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 9342cb0fcd2484ed2e5c69b2d9fdd2268d997788..a10eb417e67ab889a3312ba4de05b372adb47fd8 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -630,6 +630,10 @@ class aPCE:
# Calculate and save the score of LOOCV
score, LCerror = self.CorrectedLeaveoutError(sparsePSI, Coeffs, ModelOutput, clf_poly)
+ # Check the convergence of noise for FastARD
+ if self.RegMethod == 'FastARD' and clf_poly.alpha_ < np.finfo(np.float32).eps:
+ score = -np.inf
+
qNormScores[qidx] = score
qAllCoeffs[str(qidx+1)] = Coeffs
qAllIndices_Sparse[str(qidx+1)] = sparseBasisIndices
@@ -681,6 +685,7 @@ class aPCE:
# ------------------ Summary of results ------------------
# Select the one with the best score and save the necessary outputs
+ print("scores:",scores)
bestIdx = np.nanargmax(scores)+1
coeffs = AllCoeffs[str(bestIdx)]
PolynomialDegrees = AllIndices_Sparse[str(bestIdx)]
@@ -857,8 +862,9 @@ class aPCE:
# Transform via Principal Component Analysis
from sklearn.decomposition import PCA as sklearnPCA
+ # from sklearn.decomposition import TruncatedSVD as TruncatedSVD
n_samples, n_features = Output.shape
- covar_matrix = sklearnPCA(n_components=None)
+ covar_matrix = sklearnPCA(n_components=n_features) #None
covar_matrix.fit(Output)
var = np.cumsum(np.round(covar_matrix.explained_variance_ratio_, decimals=5)*100)
@@ -866,9 +872,9 @@ class aPCE:
selected_n_components = np.where(var>=99.999)[0][0] + 1
except:
selected_n_components = min(n_samples, n_features)
-
- nComponents = min(n_samples, n_features, selected_n_components)
+ nComponents = min(n_samples, n_features, selected_n_components)
+
pca = sklearnPCA(n_components=nComponents)
OutputMatrix = pca.fit_transform(Output)
diff --git a/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py b/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py
index 22b758590a561ef6800386636dcf63de7e96222b..1e45c64a1a853eccc394f7c554df67b212c94b41 100644
--- a/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py
+++ b/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py
@@ -226,22 +226,20 @@ if __name__ == "__main__":
BayesOpts.Name = 'Calib'
BayesOpts.emulator = True
- # MCMC solver
- BayesOpts.SamplingMethod = 'MCMC'
- BayesOpts.MCMCnSteps = 1000
-
- # Evaluation of surrogate model predictions
+ # BME Bootstrap
+ BayesOpts.Bootstrap = True
BayesOpts.NrofSamples = 100000
-
- # Bootstrap for BME calulations
- BayesOpts.Bootstrap = False
- BayesOpts.BootstrapItrNr = 10
+ BayesOpts.BootstrapItrNr = 10
BayesOpts.BootstrapNoise = 0.05
+ # Select the inference method
+ BayesOpts.SamplingMethod = 'MCMC'
+ BayesOpts.MCMCnSteps = 1000
+ BayesOpts.MultiProcessMCMC = True
+
BayesOpts.PlotPostDist = True
BayesOpts.PlotPostPred = True
-
# ----- Define the discrepancy model -------
obsData = pd.DataFrame(Model.Observations, columns=Model.Output.Names)
# (Option B)
@@ -289,20 +287,22 @@ if __name__ == "__main__":
BayesOpts_Valid = BayesInference(PCEModel)
BayesOpts_Valid.Name = 'Valid'
- BayesOpts_Valid.SamplingMethod = 'MCMC'
- BayesOpts_Valid.MCMCnSteps = 2000
- BayesOpts_Valid.Samples = Bayes_Calib.Posterior_df
BayesOpts_Valid.emulator = True
- # Bootstrap for BME calulations
- BayesOpts_Valid.Bootstrap = False
- BayesOpts_Valid.BootstrapItrNr = 20
+ # BME Bootstrap
+ BayesOpts_Valid.Bootstrap = True
+ BayesOpts_Valid.Samples = Bayes_Calib.Posterior_df
+ BayesOpts_Valid.BootstrapItrNr = 10
BayesOpts_Valid.BootstrapNoise = 0.05
+ # Select the inference method
+ BayesOpts_Valid.SamplingMethod = 'MCMC'
+ BayesOpts_Valid.MCMCnSteps = 1000
+ BayesOpts_Valid.MultiProcessMCMC = True
+
BayesOpts_Valid.PlotPostDist = True
BayesOpts_Valid.PlotPostPred = True
-
# ----- Define the discrepancy model -------
obsDataValid = pd.DataFrame(Model.ObservationsValid, columns=Model.Output.Names)
# (Option B)
diff --git a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
index b0ee08ee2dc91b3a8438fce7902ab830a26d7599..bf4ccf4c69648d521d2f9d8386413f867f57adbc 100755
--- a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
+++ b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
@@ -39,7 +39,7 @@ from BayesInference.BayesInference import BayesInference, Discrepancy
if __name__ == "__main__":
# Number of parameters
- ndim = 2 # 2, 10
+ ndim = 10 # 2, 10
#=====================================================
#============= COMPUTATIONAL MODEL ================
@@ -93,10 +93,10 @@ if __name__ == "__main__":
# error (or the highest score=1-LOO)estimator is chosen as the final
# metamodel.
MetaModelOpts.MinPceDegree = 1 #12
- MetaModelOpts.MaxPceDegree = 10 #12
+ MetaModelOpts.MaxPceDegree = 5 #12
# q-quasi-norm 0<q<1 (default=1)
- MetaModelOpts.q = 1.0 if ndim<5 else 0.65
+ MetaModelOpts.q = 1.0 if ndim<5 else 0.75
# Select the sparse least-square minimization method for
# the PCE coefficients calculation:
@@ -117,14 +117,14 @@ if __name__ == "__main__":
MetaModelOpts.addExpDesign()
# One-shot (normal) or Sequential Adaptive (sequential) Design
- MetaModelOpts.ExpDesign.Method = 'sequential'
- NrofInitSamples = 20 #75
+ MetaModelOpts.ExpDesign.Method = 'normal'
+ NrofInitSamples = 200 #75
MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
# Sampling methods
# 1) random 2) latin_hypercube 3) sobol 4) halton 5) hammersley 6) chebyshev(FT)
# 7) korobov 8) grid(FT) 9) nested_grid(FT) 10)user
- MetaModelOpts.ExpDesign.SamplingMethod = 'random'
+ MetaModelOpts.ExpDesign.SamplingMethod = 'latin_hypercube'
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.NrofNewSample = 1
@@ -221,21 +221,24 @@ if __name__ == "__main__":
# Compute and print RMSE error
PostPCE.relErrorPCEModel(nSamples=3000)
-
+ sys.exit('STOP!!')
#=====================================================
#======== Bayesian inference with Emulator ==========
#=====================================================
BayesOpts = BayesInference(PCEModel)
- # BayesOpts.SamplingMethod = 'MCMC'
BayesOpts.emulator = True
- # Bootstrap
- BayesOpts.Bootstrap = False
- BayesOpts.BootstrapItrNr = 1
+ # BME Bootstrap
+ BayesOpts.Bootstrap = True
+ BayesOpts.NrofSamples = 100000
+ BayesOpts.BootstrapItrNr = 10
BayesOpts.BootstrapNoise = 0.05
- # Evaluation of surrogate model predictions
- BayesOpts.NrofSamples = 100000
+ # Select the inference method
+ BayesOpts.SamplingMethod = 'MCMC'
+ BayesOpts.MCMCnSteps = 1000
+ BayesOpts.MultiProcessMCMC = True
+
BayesOpts.PlotPostPred = True