From 90f74934d6ea0187bf1d9dc353ecfa10e274b81c Mon Sep 17 00:00:00 2001 From: darcones Date: Fri, 19 Jul 2024 12:34:05 +0200 Subject: [PATCH] Fixed bug var to std from pce --- probeye/inference/koh/likelihood_models.py | 14 +++++++++----- probeye/inference/koh/solver.py | 4 ++-- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/probeye/inference/koh/likelihood_models.py b/probeye/inference/koh/likelihood_models.py index efccb98..cf27349 100644 --- a/probeye/inference/koh/likelihood_models.py +++ b/probeye/inference/koh/likelihood_models.py @@ -95,8 +95,11 @@ def loglike( std_model, std_meas, stds_are_scalar = self.std_values(prms) variance = np.power(std_model, 2) n = len(residual_vector) - ll = 0 + # Original likelihood + # ll = -1 / 2 * np.log(2 * np.pi * self.tolerance**2) # ll -= 0.5 / self.tolerance**2 * np.sum(np.square(residual_vector)+np.square(response_vector[1]-self.gamma*np.abs(residual_vector))) + + # Noise-corrected likelihood if std_meas is not None: variance += np.power(std_meas, 2) if stds_are_scalar: @@ -152,7 +155,7 @@ def loglike( variance += np.power(std_meas, 2) if stds_are_scalar: # ll -=0.5 * np.log(2 * np.pi / n * variance) - ll -=0.5 * np.log(2 * np.pi / n * population_variance) + ll -= 0.5 * np.log(2 * np.pi / n * population_variance) ll -= 0.5 * n / population_variance * np.square(mean_residual) # ll -= 0.5 * n / variance * np.square(mean_residual) ll -= 0.5 * n * sample_variance / population_variance @@ -237,9 +240,10 @@ def loglike( if std_meas is not None: variance += np.power(std_meas, 2) if stds_are_scalar: - ll -=0.5 * n * np.log(2 * np.pi ) - ll -= np.sum(np.log(sigma_model_sample)) - ll -= np.sum(np.square(residual_vector)/(2*np.square(sigma_model_sample))) + # ll -=0.5 * n * np.log(2 * np.pi ) + # ll -= np.sum(np.log(sigma_model_sample)) + # ll -= np.sum(np.square(residual_vector)/(2*np.square(sigma_model_sample))) + ll-= 0.5 * np.sum(np.square(np.divide(residual_vector , sigma_model_sample)) + np.log(2 * np.pi * np.square(sigma_model_sample ))) return ll diff --git a/probeye/inference/koh/solver.py b/probeye/inference/koh/solver.py index d51134a..b08ddbe 100644 --- a/probeye/inference/koh/solver.py +++ b/probeye/inference/koh/solver.py @@ -501,8 +501,8 @@ def evaluate_model_response( model_response_dict = forward_model(inp) model_response_vector, dist = vectorize_tuple_pce_dict(model_response_dict) mean_response_vector = np.array([chaospy.E(response, dist) for response in model_response_vector]).flatten() - var_response_vector = np.array([chaospy.Var(response, dist) for response in model_response_vector]).flatten() - model_response_vector = np.array([mean_response_vector, var_response_vector]) + std_response_vector = np.array([chaospy.Std(response, dist) for response in model_response_vector]).flatten() + model_response_vector = np.array([mean_response_vector, std_response_vector]) # compute the residuals by comparing to the experimental response exp_response_dict = forward_model.output_from_experiments[experiment_name]