Fixed bug var to std from pce

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
 
 

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]