Merge branch 'develop' into documentation

.github/workflows/flake8.yml

@@ -21,6 +21,6 @@ jobs:
       - name: flake8 Lint
         uses: py-actions/flake8@v1
         with:
-          ignore: "E501"
+          ignore: "E501,W605"
           exclude: "__init__.py, input/__init__.py"
           path: "pyerrors"

CONTRIBUTING.md

@@ -30,6 +30,6 @@ pytest --cov=pyerrors --cov-report html
 ```
 The linter `flake8` is executed with the command
 ```
-flake8 --ignore=E501 --exclude=__init__.py pyerrors
+flake8 --ignore=E501,W605 --exclude=__init__.py pyerrors
 ```
 Please make sure that all tests are passed for a new pull requests.

pyerrors/fits.py

@@ -239,7 +239,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs):
         if kwargs.get('covariance') is not None:
             cov = kwargs.get('covariance')
         else:
-            cov = covariance_matrix(np.concatenate((y, x.ravel())))
+            cov = covariance(np.concatenate((y, x.ravel())))
 
         number_of_x_parameters = int(m / x_f.shape[-1])
 
@@ -455,7 +455,7 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
         x0 = [0.1] * n_parms
 
     if kwargs.get('correlated_fit') is True:
-        cov = covariance_matrix(y)
+        cov = covariance(y)
         covdiag = np.diag(1. / np.sqrt(np.diag(cov)))
         corr = np.copy(cov)
         for i in range(len(y)):
@@ -527,7 +527,7 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
     if kwargs.get('expected_chisquare') is True:
         if kwargs.get('correlated_fit') is not True:
             W = np.diag(1 / np.asarray(dy_f))
-            cov = covariance_matrix(y)
+            cov = covariance(y)
             A = W @ jacobian(func)(fit_result.x, x)
             P_phi = A @ np.linalg.inv(A.T @ A) @ A.T
             expected_chisquare = np.trace((np.identity(x.shape[-1]) - P_phi) @ W @ cov @ W)
@@ -651,33 +651,9 @@ def residual_plot(x, y, func, fit_res):
     plt.draw()
 
 
-def covariance_matrix(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):
-        for j, jtem in enumerate(y[:i + 1]):
-            if i == j:
-                cov[i, j] = item.dvalue ** 2
-            else:
-                cov[i, j] = covariance(item, jtem)
-    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):
     """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta."""
-    cov = covariance_matrix(beta)
+    cov = covariance(beta)
     if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
         warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
 

pyerrors/obs.py

@@ -1332,21 +1332,60 @@ def correlate(obs_a, obs_b):
     return o
 
 
-def covariance(obs1, obs2, correlation=False, **kwargs):
+def covariance(obs, visualize=False, correlation=False, **kwargs):
-    """Calculates the covariance of two observables.
+    """Calculates the covariance matrix of a set of observables.
 
-    covariance(obs, obs) is equal to obs.dvalue ** 2
+    The gamma method has to be applied first to all observables.
-    The gamma method has to be applied first to both observables.
-
-    If abs(covariance(obs1, obs2)) > obs1.dvalue * obs2.dvalue, the covariance
-    is constrained to the maximum value.
 
     Parameters
     ----------
+    obs : list or numpy.ndarray
+        List or one dimensional array of Obs
+    visualize : bool
+        If True plots the corresponding normalized correlation matrix (default False).
     correlation : bool
-        if true the correlation instead of the covariance is returned (default False)
+        If True the correlation instead of the covariance is returned (default False).
+
+    Notes
+    -----
+    The covariance is estimated by calculating the correlation matrix assuming no autocorrelation and then rescaling the correlation matrix by the full errors including the previous gamma method estimate for the autocorrelation of the observables. For observables defined on a single ensemble this is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
+    $$
+    \tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}
+    $$
+    This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
     """
 
+    length = len(obs)
+    cov = np.zeros((length, length))
+    for i in range(length):
+        for j in range(i, length):
+            cov[i, j] = _covariance_element(obs[i], obs[j])
+    cov = cov + cov.T - np.diag(np.diag(cov))
+
+    corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+
+    errors = [o.dvalue for o in obs]
+    cov = np.diag(errors) @ corr @ np.diag(errors)
+
+    eigenvalues = np.linalg.eigh(cov)[0]
+    if not np.all(eigenvalues >= 0):
+        warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
+
+    if visualize:
+        plt.matshow(corr, vmin=-1, vmax=1)
+        plt.set_cmap('RdBu')
+        plt.colorbar()
+        plt.draw()
+
+    if correlation is True:
+        return corr
+    else:
+        return cov
+
+
+def _covariance_element(obs1, obs2):
+    """Estimates the covariance of two Obs objects, neglecting autocorrelations."""
+
     def expand_deltas(deltas, idx, shape, new_idx):
         """Expand deltas defined on idx to a contiguous range [new_idx[0], new_idx[-1]].
            New, empy entries are filled by 0. If idx and new_idx are of type range, the smallest
@@ -1369,29 +1408,18 @@ def covariance(obs1, obs2, correlation=False, **kwargs):
             ret[idx[i] - new_idx[0]] = deltas[i]
         return ret
 
-    def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx, w_max):
+    def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx):
-        gamma = np.zeros(w_max)
         deltas1 = expand_deltas(deltas1, idx1, len(idx1), new_idx)
         deltas2 = expand_deltas(deltas2, idx2, len(idx2), new_idx)
-        new_shape = len(deltas1)
+        return np.sum(deltas1 * deltas2)
-        max_gamma = min(new_shape, w_max)
-        # The padding for the fft has to be even
-        padding = new_shape + max_gamma + (new_shape + max_gamma) % 2
-        gamma[:max_gamma] += (np.fft.irfft(np.fft.rfft(deltas1, padding) * np.conjugate(np.fft.rfft(deltas2, padding)))[:max_gamma] + np.fft.irfft(np.fft.rfft(deltas2, padding) * np.conjugate(np.fft.rfft(deltas1, padding)))[:max_gamma]) / 2.0
-
-        return gamma
 
     if set(obs1.names).isdisjoint(set(obs2.names)):
-        return 0.
+        return 0.0
 
     if not hasattr(obs1, 'e_dvalue') or not hasattr(obs2, 'e_dvalue'):
         raise Exception('The gamma method has to be applied to both Obs first.')
 
-    dvalue = 0
+    dvalue = 0.0
-    e_gamma = {}
-    e_dvalue = {}
-    e_n_tauint = {}
-    e_rho = {}
 
     for e_name in obs1.mc_names:
 
@@ -1399,52 +1427,32 @@ def covariance(obs1, obs2, correlation=False, **kwargs):
             continue
 
         idl_d = {}
-        r_length = []
         for r_name in obs1.e_content[e_name]:
             if r_name not in obs2.e_content[e_name]:
                 continue
             idl_d[r_name] = _merge_idx([obs1.idl[r_name], obs2.idl[r_name]])
-            if isinstance(idl_d[r_name], range):
-                r_length.append(len(idl_d[r_name]))
-            else:
-                r_length.append((idl_d[r_name][-1] - idl_d[r_name][0] + 1))
 
-        if not r_length:
+        gamma = 0.0
-            return 0.
-
-        w_max = max(r_length) // 2
-        e_gamma[e_name] = np.zeros(w_max)
 
         for r_name in obs1.e_content[e_name]:
             if r_name not in obs2.e_content[e_name]:
                 continue
-            e_gamma[e_name] += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name], w_max)
+            gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name])
 
-        if np.all(e_gamma[e_name] == 0.0):
+        if gamma == 0.0:
             continue
 
-        e_shapes = []
+        gamma_div = 0.0
-        for r_name in obs1.e_content[e_name]:
-            e_shapes.append(obs1.shape[r_name])
-        gamma_div = np.zeros(w_max)
         e_N = 0
         for r_name in obs1.e_content[e_name]:
             if r_name not in obs2.e_content[e_name]:
                 continue
-            gamma_div += calc_gamma(np.ones(obs1.shape[r_name]), np.ones(obs2.shape[r_name]), obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name], w_max)
+            gamma_div += calc_gamma(np.ones(obs1.shape[r_name]), np.ones(obs2.shape[r_name]), obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name])
             e_N += np.sum(np.ones_like(idl_d[r_name]))
-        e_gamma[e_name] /= gamma_div[:w_max]
+        gamma /= gamma_div
 
-        e_rho[e_name] = e_gamma[e_name][:w_max] / e_gamma[e_name][0]
-        e_n_tauint[e_name] = np.cumsum(np.concatenate(([0.5], e_rho[e_name][1:])))
-        # Make sure no entry of tauint is smaller than 0.5
-        e_n_tauint[e_name][e_n_tauint[e_name] < 0.5] = 0.500000000001
-
-        window = min(obs1.e_windowsize[e_name], obs2.e_windowsize[e_name])
         # Bias correction hep-lat/0306017 eq. (49)
-        e_dvalue[e_name] = 2 * (e_n_tauint[e_name][window] + obs1.tau_exp[e_name] * np.abs(e_rho[e_name][window + 1])) * (1 + (2 * window + 1) / e_N) * e_gamma[e_name][0] / e_N
+        dvalue += (1 + 1 / e_N) * gamma / e_N
-
-        dvalue += e_dvalue[e_name]
 
     for e_name in obs1.cov_names:
 
@@ -1453,12 +1461,6 @@ def covariance(obs1, obs2, correlation=False, **kwargs):
 
         dvalue += float(np.dot(np.transpose(obs1.covobs[e_name].grad), np.dot(obs1.covobs[e_name].cov, obs2.covobs[e_name].grad)))
 
-    if np.abs(dvalue / obs1.dvalue / obs2.dvalue) > 1.0:
-        dvalue = np.sign(dvalue) * obs1.dvalue * obs2.dvalue
-
-    if correlation:
-        dvalue = dvalue / obs1.dvalue / obs2.dvalue
-
     return dvalue
 
 

tests/covobs_test.py

@@ -39,8 +39,8 @@ def test_covobs():
         assert(np.isclose(oc.value, op.value, rtol=1e-14, atol=1e-14))
 
     [o.gamma_method() for o in cl]
-    assert(pe.covariance(cl[0], cl[1]) == cov[0][1])
+    assert(np.isclose(pe.covariance([cl[0], cl[1]])[0, 1], cov[0][1]))
-    assert(pe.covariance(cl[0], cl[1]) == cov[1][0])
+    assert(np.isclose(pe.covariance([cl[0], cl[1]])[0, 1], cov[1][0]))
 
     do = cl[0] * cl[1]
     assert(np.array_equal(do.covobs['rAP'].grad, np.transpose([pi[1], pi[0]]).reshape(2, 1)))
@@ -82,6 +82,19 @@ def test_covobs_init():
     covobs = pe.cov_Obs([1, 2], np.array([[0.21, 0.2], [0.2, 0.21]]), 'test')
 
 
+def test_covobs_covariance():
+    a = pe.cov_Obs(2.47, 0.03 ** 2, "Cov_obs 1")
+    b = pe.cov_Obs(-4.3, 0.335 ** 2, "Cov_obs 2")
+
+    x = [a + b, a - b]
+    [o.gamma_method() for o in x]
+
+    covariance = pe.covariance(x)
+
+    assert np.isclose(covariance[0, 0], covariance[1, 1])
+    assert np.isclose(covariance[0, 1], a.dvalue ** 2 - b.dvalue ** 2)
+
+
 def test_covobs_exceptions():
     with pytest.raises(Exception):
         covobs = pe.cov_Obs(0.1, [[0.1, 0.2], [0.1, 0.2]], 'test')

tests/fits_test.py

@@ -60,7 +60,7 @@ def test_least_squares():
         beta[i].gamma_method(S=1.0)
         assert math.isclose(beta[i].value, popt[i], abs_tol=1e-5)
         assert math.isclose(pcov[i, i], beta[i].dvalue ** 2, abs_tol=1e-3)
-    assert math.isclose(pe.covariance(beta[0], beta[1]), pcov[0, 1], abs_tol=1e-3)
+    assert math.isclose(pe.covariance([beta[0], beta[1]])[0, 1], pcov[0, 1], abs_tol=1e-3)
 
     chi2_pyerrors = np.sum(((f(x, *[o.value for o in beta]) - y) / yerr) ** 2) / (len(x) - 2)
     chi2_scipy = np.sum(((f(x, *popt) - y) / yerr) ** 2) / (len(x) - 2)
@@ -81,7 +81,7 @@ def test_least_squares():
         betac[i].gamma_method(S=1.0)
         assert math.isclose(betac[i].value, popt[i], abs_tol=1e-5)
         assert math.isclose(pcov[i, i], betac[i].dvalue ** 2, abs_tol=1e-3)
-    assert math.isclose(pe.covariance(betac[0], betac[1]), pcov[0, 1], abs_tol=1e-3)
+    assert math.isclose(pe.covariance([betac[0], betac[1]])[0, 1], pcov[0, 1], abs_tol=1e-3)
 
 
 def test_alternative_solvers():
@@ -195,7 +195,7 @@ def test_total_least_squares():
         beta[i].gamma_method(S=1.0)
         assert math.isclose(beta[i].value, output.beta[i], rel_tol=1e-5)
         assert math.isclose(output.cov_beta[i, i], beta[i].dvalue ** 2, rel_tol=2.5e-1), str(output.cov_beta[i, i]) + ' ' + str(beta[i].dvalue ** 2)
-    assert math.isclose(pe.covariance(beta[0], beta[1]), output.cov_beta[0, 1], rel_tol=2.5e-1)
+    assert math.isclose(pe.covariance([beta[0], beta[1]])[0, 1], output.cov_beta[0, 1], rel_tol=2.5e-1)
 
     out = pe.total_least_squares(ox, oy, func, const_par=[beta[1]])
 
@@ -218,7 +218,7 @@ def test_total_least_squares():
         betac[i].gamma_method(S=1.0)
         assert math.isclose(betac[i].value, output.beta[i], rel_tol=1e-3)
         assert math.isclose(output.cov_beta[i, i], betac[i].dvalue ** 2, rel_tol=2.5e-1), str(output.cov_beta[i, i]) + ' ' + str(betac[i].dvalue ** 2)
-    assert math.isclose(pe.covariance(betac[0], betac[1]), output.cov_beta[0, 1], rel_tol=2.5e-1)
+    assert math.isclose(pe.covariance([betac[0], betac[1]])[0, 1], output.cov_beta[0, 1], rel_tol=2.5e-1)
 
     outc = pe.total_least_squares(oxc, oyc, func, const_par=[betac[1]])
 
@@ -233,7 +233,7 @@ def test_total_least_squares():
         betac[i].gamma_method(S=1.0)
         assert math.isclose(betac[i].value, output.beta[i], rel_tol=1e-3)
         assert math.isclose(output.cov_beta[i, i], betac[i].dvalue ** 2, rel_tol=2.5e-1), str(output.cov_beta[i, i]) + ' ' + str(betac[i].dvalue ** 2)
-    assert math.isclose(pe.covariance(betac[0], betac[1]), output.cov_beta[0, 1], rel_tol=2.5e-1)
+    assert math.isclose(pe.covariance([betac[0], betac[1]])[0, 1], output.cov_beta[0, 1], rel_tol=2.5e-1)
 
     outc = pe.total_least_squares(oxc, oy, func, const_par=[betac[1]])
 

tests/obs_test.py

@@ -246,17 +246,6 @@ def test_gamma_method_kwargs():
     assert my_obs.N_sigma['ens'] == pe.Obs.N_sigma_global
 
 
-def test_covariance_is_variance():
-    value = np.random.normal(5, 10)
-    dvalue = np.abs(np.random.normal(0, 1))
-    test_obs = pe.pseudo_Obs(value, dvalue, 't')
-    test_obs.gamma_method()
-    assert np.abs(test_obs.dvalue ** 2 - pe.covariance(test_obs, test_obs)) <= 10 * np.finfo(np.float64).eps
-    test_obs = test_obs + pe.pseudo_Obs(value, dvalue, 'q', 200)
-    test_obs.gamma_method()
-    assert np.abs(test_obs.dvalue ** 2 - pe.covariance(test_obs, test_obs)) <= 10 * np.finfo(np.float64).eps
-
-
 def test_fft():
     value = np.random.normal(5, 100)
     dvalue = np.abs(np.random.normal(0, 5))
@@ -268,21 +257,6 @@ def test_fft():
     assert np.abs(test_obs1.dvalue - test_obs2.dvalue) <= 10 * max(test_obs1.dvalue, test_obs2.dvalue) * np.finfo(np.float64).eps
 
 
-def test_covariance_symmetry():
-    value1 = np.random.normal(5, 10)
-    dvalue1 = np.abs(np.random.normal(0, 1))
-    test_obs1 = pe.pseudo_Obs(value1, dvalue1, 't')
-    test_obs1.gamma_method()
-    value2 = np.random.normal(5, 10)
-    dvalue2 = np.abs(np.random.normal(0, 1))
-    test_obs2 = pe.pseudo_Obs(value2, dvalue2, 't')
-    test_obs2.gamma_method()
-    cov_ab = pe.covariance(test_obs1, test_obs2)
-    cov_ba = pe.covariance(test_obs2, test_obs1)
-    assert np.abs(cov_ab - cov_ba) <= 10 * np.finfo(np.float64).eps
-    assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float64).eps)
-
-
 def test_gamma_method_uncorrelated():
     # Construct pseudo Obs with random shape
     value = np.random.normal(5, 10)
@@ -620,6 +594,18 @@ def test_gamma_method_irregular():
     ol = [a, b]
     o = (ol[0] - ol[1]) / (ol[1])
 
+
+def test_covariance_is_variance():
+    value = np.random.normal(5, 10)
+    dvalue = np.abs(np.random.normal(0, 1))
+    test_obs = pe.pseudo_Obs(value, dvalue, 't')
+    test_obs.gamma_method()
+    assert np.isclose(test_obs.dvalue ** 2, pe.covariance([test_obs, test_obs])[0, 1])
+    test_obs = test_obs + pe.pseudo_Obs(value, dvalue, 'q', 200)
+    test_obs.gamma_method()
+    assert np.isclose(test_obs.dvalue ** 2, pe.covariance([test_obs, test_obs])[0, 1])
+
+
 def test_covariance_symmetry():
     value1 = np.random.normal(5, 10)
     dvalue1 = np.abs(np.random.normal(0, 1))
@@ -629,9 +615,9 @@ def test_covariance_symmetry():
     dvalue2 = np.abs(np.random.normal(0, 1))
     test_obs2 = pe.pseudo_Obs(value2, dvalue2, 't')
     test_obs2.gamma_method()
-    cov_ab = pe.covariance(test_obs1, test_obs2)
+    cov_ab = pe.covariance([test_obs1, test_obs2])[0, 1]
-    cov_ba = pe.covariance(test_obs2, test_obs1)
+    cov_ba = pe.covariance([test_obs2, test_obs1])[0, 1]
-    assert np.abs(cov_ab - cov_ba) <= 10 * np.finfo(np.float64).eps
+    assert np.isclose(cov_ab, cov_ba)
     assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float64).eps)
 
     N = 100
@@ -643,14 +629,73 @@ def test_covariance_symmetry():
     idx = [i + 1 for i in range(len(configs)) if configs[i] == 1]
     a = pe.Obs([zero_arr], ['t'], idl=[idx])
     a.gamma_method()
-    assert np.isclose(a.dvalue**2, pe.covariance(a, a), atol=100, rtol=1e-4)
+    assert np.isclose(a.dvalue ** 2, pe.covariance([a, a])[0, 1], atol=100, rtol=1e-4)
 
-    cov_ab = pe.covariance(test_obs1, a)
+    cov_ab = pe.covariance([test_obs1, a])[0, 1]
-    cov_ba = pe.covariance(a, test_obs1)
+    cov_ba = pe.covariance([a, test_obs1])[0, 1]
     assert np.abs(cov_ab - cov_ba) <= 10 * np.finfo(np.float64).eps
     assert np.abs(cov_ab) < test_obs1.dvalue * a.dvalue * (1 + 10 * np.finfo(np.float64).eps)
 
 
+def test_covariance_sum():
+    length = 2
+    t_fac = 0.4
+    tt = pe.misc.gen_correlated_data(np.zeros(length), 0.99 * np.ones((length, length)) + 0.01 * np.diag(np.ones(length)), 'test', tau=0.5 + t_fac * np.random.rand(length), samples=1000)
+    [o.gamma_method(S=0) for o in tt]
+
+    t_cov = pe.covariance(tt)
+
+    my_sum = tt[0] + tt[1]
+    my_sum.gamma_method(S=0)
+    e_cov = (my_sum.dvalue ** 2 - tt[0].dvalue ** 2 - tt[1].dvalue ** 2) / 2
+
+    assert np.isclose(e_cov, t_cov[0, 1])
+
+
+def test_covariance_positive_semidefinite():
+    length = 64
+    t_fac = 1.5
+    tt = pe.misc.gen_correlated_data(np.zeros(length), 0.99999 * np.ones((length, length)) + 0.00001 * np.diag(np.ones(length)), 'test', tau=0.5 + t_fac * np.random.rand(length), samples=1000)
+    [o.gamma_method() for o in tt]
+    cov = pe.covariance(tt)
+    assert np.all(np.linalg.eigh(cov)[0] >= -1e-15)
+
+
+def test_covariance_factorizing():
+    length = 2
+    t_fac = 1.5
+
+    tt = pe.misc.gen_correlated_data(np.zeros(length), 0.75 * np.ones((length, length)) + 0.8 * np.diag(np.ones(length)), 'test', tau=0.5 + t_fac * np.random.rand(length), samples=1000)
+    [o.gamma_method() for o in tt]
+
+    mt0 = -tt[0]
+    mt0.gamma_method()
+
+    assert np.isclose(pe.covariance([mt0, tt[1]])[0, 1], -pe.covariance(tt)[0, 1])
+
+
+def test_covariance_alternation():
+    length = 12
+    t_fac = 2.5
+
+    tt1 = pe.misc.gen_correlated_data(np.zeros(length), -0.00001 * np.ones((length, length)) + 0.002 * np.diag(np.ones(length)), 'test', tau=0.5 + t_fac * np.random.rand(length), samples=88)
+    tt2 = pe.misc.gen_correlated_data(np.zeros(length), 0.9999 * np.ones((length, length)) + 0.0001 * np.diag(np.ones(length)), 'another_test|r0', tau=0.7 + t_fac * np.random.rand(length), samples=73)
+    tt3 = pe.misc.gen_correlated_data(np.zeros(length), 0.9999 * np.ones((length, length)) + 0.0001 * np.diag(np.ones(length)), 'another_test|r1', tau=0.7 + t_fac * np.random.rand(length), samples=91)
+
+    tt = np.array(tt1) + (np.array(tt2) + np.array(tt3))
+    tt *= np.resize([1, -1], length)
+
+    [o.gamma_method() for o in tt]
+    cov = pe.covariance(tt, True)
+
+    assert np.all(np.linalg.eigh(cov)[0] > -1e-15)
+
+
+def test_covariance_correlation():
+    test_obs = pe.pseudo_Obs(-4, 0.8, 'test', samples=784)
+    assert np.allclose(pe.covariance([test_obs, test_obs, test_obs], correlation=True), np.ones((3, 3)))
+
+
 def test_empty_obs():
     o = pe.Obs([np.random.rand(100)], ['test'])
     q = o + pe.Obs([], [])