feat: a warning is now issued when an estimated covariance matrix is not

positive semi-definite. Docstrings extended.

pyerrors/fits.py

@@ -652,7 +652,13 @@ def residual_plot(x, y, func, fit_res):
 
 
 def covariance_matrix(y):
-    """Returns the covariance matrix of y."""
+    """Returns the covariance matrix of y.
+
+    Parameters
+    ----------
+    y : list or numpy.ndarray
+        List or one dimensional array of Obs
+    """
     length = len(y)
     cov = np.zeros((length, length))
     for i, item in enumerate(y):
@@ -661,7 +667,12 @@ def covariance_matrix(y):
                 cov[i, j] = item.dvalue ** 2
             else:
                 cov[i, j] = covariance(item, jtem)
-    return cov + cov.T - np.diag(np.diag(cov))
+    cov = cov + cov.T - np.diag(np.diag(cov))
+    eigenvalues = np.linalg.eigh(cov)[0]
+    if not np.all(eigenvalues >= 0):
+        warnings.warn("Covariance matrix is not positive semi-definite", RuntimeWarning)
+        print("Eigenvalues of the covariance matrix:", eigenvalues)
+    return cov
 
 
 def error_band(x, func, beta):

pyerrors/obs.py

@@ -1341,10 +1341,10 @@ def covariance(obs1, obs2, correlation=False, **kwargs):
     If abs(covariance(obs1, obs2)) > obs1.dvalue * obs2.dvalue, the covariance
     is constrained to the maximum value.
 
-    Keyword arguments
-    -----------------
-    correlation -- if true the correlation instead of the covariance is
-                   returned (default False)
+    Parameters
+    ----------
+    correlation : bool
+        if true the correlation instead of the covariance is returned (default False)
     """
 
     def expand_deltas(deltas, idx, shape, new_idx):