fix pva and add a seed for the tests (#709)

smt/utils/misc.py

@@ -120,31 +120,84 @@ def compute_relative_error(sm, xe=None, ye=None, kx=None):
 
 
 def compute_pva(sm, xe, ye):
-    ye = ye.reshape((xe.shape[0], 1))
-    N = len(ye)
-    ye2 = sm.predict_values(xe)
-    variance = sm.predict_variances(xe)
+    """
+    Compute the Predictive Variance Adequacy (PVA) for a surrogate model.
+
+    Parameters:
+    - sm: The surrogate model object, expected to have `predict_values` and `predict_variances` methods.
+    - xe: Input data for evaluation (N x d array).
+    - ye: True output values (N x 1 array or equivalent).
+
+    Returns:
+    - pva: Predictive Variance Adequacy score (float).
+    """
+    ye = ye.reshape((xe.shape[0], 1))  # Ensure `ye` is column vector
+    Nb = len(ye)  # Number of data points
+
+    # Predicted values and variances
+    ye_pred = sm.predict_values(xe)  # Predicted values (N x 1 array)
+    variance = sm.predict_variances(xe)  # Predicted variances (N x 1 array)
+
+    # Calculate squared error normalized by variance
+    error = ((ye_pred - ye) ** 2) / variance
+
+    # Compute PVA with logarithm
+    pva = np.abs(np.log(np.sum(error) / Nb))
 
-    error = (ye2 - ye) ** 2 / variance
-    pva = np.sum(error) / N
     return pva
 
 
 def compute_rmse(sm, xe, ye):
-    ye = ye.reshape((xe.shape[0], 1))
-    N = len(ye)
-    ye2 = sm.predict_values(xe)
-    rmse = np.sqrt(np.sum((ye2 - ye) ** 2) / N)
+    """
+    Compute the Root Mean Square Error (RMSE) for a surrogate model.
+
+    Parameters:
+    - sm: The surrogate model object, expected to have a `predict_values` method.
+    - xe: Input data for evaluation (N x d array).
+    - ye: True output values (N x 1 array or equivalent).
+
+    Returns:
+    - rmse: Root Mean Square Error (float).
+    """
+    ye = ye.reshape((xe.shape[0], 1))  # Ensure `ye` is a column vector
+
+    # Predicted values
+    ye_pred = sm.predict_values(xe)  # Predicted values (N x 1 array)
+
+    # Compute RMSE
+    mse = np.mean((ye_pred - ye) ** 2)  # Mean Squared Error
+    rmse = np.sqrt(mse)  # Root Mean Squared Error
+
     return rmse
 
 
 def compute_q2(sm, xe, ye):
-    ye = ye.reshape((xe.shape[0], 1))
-    N = len(ye)
-    square_rmse = compute_rmse(sm, xe, ye) ** 2
+    """
+    Compute the Q^2 validation criterion for a surrogate model.
+
+    Parameters:
+    - sm: The surrogate model object, expected to have a `predict_values` method.
+    - xe: Input data for evaluation (N x d array).
+    - ye: True output values (N x 1 array or equivalent).
+
+    Returns:
+    - Q2: Predictive coefficient of determination (float).
+    """
+    ye = ye.reshape((xe.shape[0], 1))  # Ensure `ye` is a column vector
+
+    # Predicted values
+    ye_pred = sm.predict_values(xe)  # Predicted values (N x 1 array)
+
+    # Mean of true output values
     ye_mean = np.mean(ye)
-    variance = np.sum((ye - ye_mean) ** 2) / N
-    Q2 = 1 - (square_rmse / variance)
+
+    # Residual Sum of Squares (RSS) and Total Sum of Squares (TSS)
+    rss = np.sum((ye - ye_pred) ** 2)  # Residual sum of squares
+    tss = np.sum((ye - ye_mean) ** 2)  # Total sum of squares
+
+    # Compute Q^2
+    Q2 = 1 - (rss / tss)
+
     return Q2
 
 

smt/utils/test/test_misc_utils.py

@@ -40,37 +40,39 @@ def prepare_tests_errors(self):
         ndim = 2
         fun = Sphere(ndim=ndim)
 
-        sampling = LHS(xlimits=fun.xlimits, criterion="m")
-        xe = sampling(20)
+        sampling = LHS(xlimits=fun.xlimits, criterion="ese", random_state=42)
+        xt = sampling(20)
+        yt = fun(xt)
+        xe = sampling(100)
         ye = fun(xe)
-        return xe, ye
+        return xt, yt, xe, ye
 
     def test_pva_error(self):
-        xe, ye = self.prepare_tests_errors()
-        sm = KRG(print_global=False)
-        sm.set_training_values(xe, ye)
+        xt, yt, xe, ye = self.prepare_tests_errors()
+        sm = KRG(print_global=False, n_start=25, random_state=42)
+        sm.set_training_values(xt, yt)
         sm.train()
 
         pva = compute_pva(sm, xe, ye)
-        self.assertLess(pva, 0.7)
+        self.assertAlmostEqual(pva, 0.1, delta=2e-2)
 
     def test_rmse_error(self):
-        xe, ye = self.prepare_tests_errors()
-        sm = KRG(print_global=False)
-        sm.set_training_values(xe, ye)
+        xt, yt, xe, ye = self.prepare_tests_errors()
+        sm = KRG(print_global=False, random_state=42)
+        sm.set_training_values(xt, yt)
         sm.train()
 
         rmse = compute_rmse(sm, xe, ye)
-        self.assertLess(rmse, 0.1)
+        self.assertAlmostEqual(rmse, 0.0, delta=1e-2)
 
     def test_q2_error(self):
-        xe, ye = self.prepare_tests_errors()
-        sm = KRG(print_global=False)
-        sm.set_training_values(xe, ye)
+        xt, yt, xe, ye = self.prepare_tests_errors()
+        sm = KRG(print_global=False, random_state=42)
+        sm.set_training_values(xt, yt)
         sm.train()
 
         q2 = compute_q2(sm, xe, ye)
-        self.assertAlmostEqual(q2, 1.0, delta=1e-3)
+        self.assertAlmostEqual(q2, 1.0, delta=1e-2)
 
 
 if __name__ == "__main__":