diff --git a/src/bayesvalidrox/bayes_inference/bayes_inference.py b/src/bayesvalidrox/bayes_inference/bayes_inference.py
index 0ecbd33b774531077e8073725ac96eac43e44a16..c8073ba76ebeb3249646e17b7f395f8a2472e2f6 100644
--- a/src/bayesvalidrox/bayes_inference/bayes_inference.py
+++ b/src/bayesvalidrox/bayes_inference/bayes_inference.py
@@ -547,16 +547,6 @@ class BayesInference:
# -------- Plot log_BME dist --------
if self.bootstrap:
- if hasattr(self, 'sigma2s'):
- sigma2 = self.sigma2s
- else:
- sigma2 = None
- # Model validation via binary hypothesis testing
- valid_outputs = self.model_validation(
- self.measured_data,
- self.Discrepancy.total_sigma2,
- sigma2=sigma2,
- std=surrError)
# Computing the TOM performance
self.log_BME_tom = stats.chi2.rvs(
@@ -652,170 +642,6 @@ class BayesInference:
return final_data
- # -------------------------------------------------------------------------
- def model_validation(self, obs_data, total_sigma2s, sigma2=None, std=None):
- """
- A Bayesian risk-based model validation based on the following paper:
- Jiang, Xiaomo, and Sankaran Mahadevan. "Bayesian validation
- assessment of multivariate computational models." Journal of
- Applied Statistics 35, no. 1 (2008): 49-65.
-
- Parameters
- ----------
- obs_data : dict
- A dictionary/dataframe containing the observation data.
- total_sigma2s : dict
- A dictionary with known values of the covariance diagonal entries,
- a.k.a sigma^2.
- sigma2 : array, optional
- An array of the sigma^2 samples, when the covariance diagonal
- entries are unknown and are being jointly inferred. The default is
- None.
- std : dict, optional
- A dictionary containing the root mean squared error as array of
- shape (n_samples, n_measurement) for each model output. The default
- is None.
-
- Returns
- -------
- final_output : dict
- It contains the following values for each output model:
-
- 1. Bayes Factor
- 2. Confidence in model acceptance (in %)
- 3. Probability of model acceptance (in %)
- """
-
- # Set some variables
- MetaModel = self.MetaModel
- Model = MetaModel.ModelObj
- output_names = Model.Output.names
-
- n_samples = 10000
- samples = MetaModel.ExpDesign.generate_samples(n_samples, 'random')
-
- # Generate model prediction
- if self.emulator:
- y_hat, y_std = MetaModel.eval_metamodel(samples=samples)
- else:
- y_hat = self._eval_model(
- samples=samples, key='Validation'
- )
-
- # Compute Bayes factor in the logarichmic scale
- final_output = {}
- for out_idx, output in enumerate(output_names):
- output_bf = {}
- ln_BF = np.zeros((self.n_bootstrap_itrs))
- # Generate data using perturb data
- perturb_data = self._perturb_data(obs_data, [output])
- n_out = y_hat[output].shape[1]
- for i, data in enumerate(perturb_data):
- # Compute \bar{Y}
- y_bar = np.mean(data - y_hat[output], axis=0)
- # Prepare \Sigma
- non_nan_indices = ~np.isnan(total_sigma2s[output])
- tot_sigma2s = total_sigma2s[output][non_nan_indices]
-
- # Add the std of the PCE is chosen as emulator.
- if self.emulator:
- if std is not None:
- tot_sigma2s += std[output]**2
-
- # Covariance Matrix
- covMatrix = np.diag(tot_sigma2s)
- covMatrix += np.cov((data - y_hat[output]).T)
-
- # Infer equal sigma2s
- try:
- sigma_2 = np.mean(sigma2[:, out_idx])
- except TypeError:
- sigma_2 = 0.0
-
- covMatrix += sigma_2 * np.eye(n_out)
-
- # Select the data points to compare
- try:
- indices = self.selected_indices[output]
- except:
- indices = list(range(n_out))
- covMatrix = np.diag(covMatrix[indices, indices])
- y_bar = y_bar[indices]
-
- # Bayes factor Eq. 11
- ln_BF[i] = -(n_samples**2/(2*n_samples+2))
- ln_BF[i] *= np.dot(np.dot(y_bar, np.linalg.pinv(covMatrix)),
- np.transpose(y_bar))
- ln_BF[i] += np.log(n_samples+1)/2
-
- output_bf['Bayes Factor'] = ln_BF
-
- # Compute confidence in model acceptance (in %) Eq. 16
- pr_H0_Y = np.exp(output_bf['Bayes Factor']) * 100
- pr_H0_Y /= (1+np.exp(output_bf['Bayes Factor']))
- output_bf['Confidence in model acceptance'] = pr_H0_Y
-
- # Compute the probability of accepting the model (γ) Eq. 19
- threshold = np.log(1) # Eq. 4
- gamma = np.sum(output_bf['Bayes Factor'] > threshold) * 100
- gamma /= self.n_bootstrap_itrs
- output_bf['Prob. of model acceptance'] = gamma
-
- # Store computed quantities for each model output
- final_output[output] = output_bf
-
- # Save the results
- with open('model_validation_output.pkl', 'wb') as output:
- joblib.dump(final_output, output, 2)
-
- # Plot the ln Bayes factors
- # Create Output directory, if it doesn't exist already.
- out_dir = f'Outputs_Bayes_{Model.name}_{self.name}'
- os.makedirs(out_dir, exist_ok=True)
-
- for out_idx, output in enumerate(output_names):
- # Create figure
- fig, ax = plt.subplots()
- # Plot the distirbution of BF
- sns.kdeplot(
- final_output[output]['Bayes Factor'],
- ax=ax, color="green", shade=True
- )
- # Plot the vertical line
- decision_threshold = 1
- plt.axvline(np.log(decision_threshold), c='red',
- ls='--')
-
- ax.set_xlabel('ln(BF)')
- ax.set_ylabel('Probability density')
-
- legend_elements = [
- Patch(facecolor='green', edgecolor='green',
- label=f'{Model.name} - {output}')
- ]
- ax.legend(handles=legend_elements)
-
- # Annotate the Prob. of model acceptance
- prob_accept = final_output[output]['Prob. of model acceptance']
- text = f"Prob. accpetance={prob_accept:.2f} %"
-
- plt.text(0.15, 0.75, text, bbox=dict(facecolor='wheat',
- edgecolor='black',
- boxstyle='round,pad=1'),
- transform=ax.transAxes)
-
- if self.emulator:
- plotname = f'/BME_hist_{Model.name}_emulator_{output}'
- else:
- plotname = f'/BME_hist_{Model.name}_{output}'
-
- # Save figure
- plt.savefig(f'./{out_dir}{plotname}.pdf', bbox_inches='tight')
- plt.show()
- plt.close()
-
- return final_output
-
# -------------------------------------------------------------------------
def _logpdf(self, x, mean, cov):
"""