From 310ca20474bcd39223db3a1d4b08eda51e641023 Mon Sep 17 00:00:00 2001
From: farid <farid.mohammadi@iws.uni-stuttgart.de>
Date: Fri, 7 Jan 2022 12:25:32 +0100
Subject: [PATCH] [Bayesinference] refactor likelihhod function.
---
.../BayesInference/BayesInference.py | 101 +++++++++---------
1 file changed, 49 insertions(+), 52 deletions(-)
diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index 1072faab2..43f4d789a 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -214,72 +214,69 @@ class BayesInference:
else:
OutputType = list(self.ReqOutputType)
- # Loop over the outputs and calculate logLikelihood
+ # Loop over the outputs
for idx, out in enumerate(OutputType):
# (Meta)Model Output
- TotalOutputs = Outputs[out]
- nout = TotalOutputs.shape[1]
+ nsamples, nout = Outputs[out].shape
- # Remove NaN
+ # Prepare data and remove NaN
try:
data = Data[out].to_numpy()[~np.isnan(Data[out])]
except:
data = Data[out][~np.isnan(Data[out])]
totalSigma2s = TotalSigma2s[out][~np.isnan(TotalSigma2s[out])][:nout]
- # Covariance Matrix
- covMatrix = totalSigma2s
-
- if Sigma2 is not None:
- # Check the type error term
- if hasattr(self, 'BiasInputs') and not hasattr(self, 'errorModel'):
- # TODO: Infer a Bias model usig Gaussian Process Regression
- BiasInputs = np.hstack((self.BiasInputs[out],TotalOutputs.T))
- # EDY = self.PCEModel.ExpDesign.Y[out]
- # BiasInputs = np.hstack((self.BiasInputs[out],EDY.T))
- sigma2 = Sigma2[:3] if idx==0 else Sigma2[3:]
- biasSigma2s = self.RBF_kernel(BiasInputs, sigma2)
- # biasSigma2s = self.RBF_kernel(self.BiasInputs[out], Sigma2)
- else:
- # Infer equal sigma2s
- try:
- sigma2 = Sigma2[:,idx]
- except:
- sigma2 = 0.0
+ # Loop over each run/sample and calculate logLikelihood
+ logliks = np.zeros(nsamples)
+ for s_idx in range(nsamples):
+
+ # Simulation run
+ TotalOutputs = Outputs[out]
+
+ # Covariance Matrix
+ covMatrix = np.diag(totalSigma2s)
- # Convert biasSigma2s to a covMatrix
- # biasSigma2s = sigma2 * np.eye(nout)
- biasSigma2s = np.multiply( sigma2.reshape(-1,1) , totalSigma2s.reshape(1,-1))
-
- # Add Sigma2 to covMatrix
- # covMatrix += biasSigma2s
- covMatrix = biasSigma2s
+ if Sigma2 is not None:
+ # Check the type error term
+ if hasattr(self, 'BiasInputs') and not hasattr(self, 'errorModel'):
+ # TODO: Infer a Bias model usig Gaussian Process Regression
+ BiasInputs = np.hstack((self.BiasInputs[out],
+ TotalOutputs[s_idx].reshape(-1,1)))
+ # EDY = self.PCEModel.ExpDesign.Y[out]
+ # BiasInputs = np.hstack((self.BiasInputs[out],EDY.T))
+ params = Sigma2[s_idx,:3] if idx==0 else Sigma2[s_idx,3:]
+ covMatrix = self.RBF_kernel(BiasInputs, params)
+ # covMatrix = self.RBF_kernel(self.BiasInputs[out], Sigma2)
+ else:
+ # Infer equal sigma2s
+ try:
+ sigma2 = Sigma2[s_idx,idx]
+ except:
+ sigma2 = 0.0
- # Add the std of the PCE is chosen as emulator.
- if self.emulator:
- stdPCE = std[out] if std is not None else np.mean(self._stdPCEPriorPred[out], axis=0)
- # Expected value of variance (Assump: i.i.d stds)
- covMatrix += stdPCE**2
+ # covMatrix = np.diag(sigma2 * np.eye(nout))
+ covMatrix = np.diag(sigma2 * totalSigma2s)
+
+ # Add the std of the PCE is chosen as emulator.
+ if self.emulator:
+ stdPCE = std[out] if std is not None else np.mean(self._stdPCEPriorPred[out], axis=0)
+ # Expected value of variance (Assump: i.i.d stds)
+ covMatrix += np.diag(stdPCE**2)
- # print(formatting_function(covMatrix))
- # Select the data points to compare
- if self.selectedIndices is not None:
- TotalOutputs = TotalOutputs[:,self.selectedIndices[out]]
- data = data[self.selectedIndices[out]]
- covMatrix = covMatrix[:,self.selectedIndices[out]]
+ # print(formatting_function(covMatrix))
+ # Select the data points to compare
+ if self.selectedIndices is not None:
+ indices = self.selectedIndices[out]
+ covMatrix = np.diag(covMatrix[indices,indices])
+ else:
+ indices = list(range(nout))
- # Compute loglikelihood
- # if TotalOutputs.shape[0] == 1:
- # logLik += self._logpdf(TotalOutputs[0], data, np.diag(covMatrix))
- if covMatrix.ndim == 1:
- logLik += stats.multivariate_normal.logpdf(TotalOutputs, data, np.diag(covMatrix))
- else:
- logliks = np.zeros(TotalOutputs.shape[0])
- for i,cov in enumerate(covMatrix):
- logliks[i] = stats.multivariate_normal.logpdf(TotalOutputs[i], data, np.diag(cov))
-
- logLik += logliks
+ # Compute loglikelihood
+ logliks[s_idx] = self._logpdf(TotalOutputs[s_idx,indices],
+ data[indices], covMatrix)
+
+ logLik += logliks
return logLik
--
GitLab