mirror of
https://github.com/fjosw/pyerrors.git
synced 2025-03-15 14:50:25 +01:00
feat: the correlation for two observables with different idls is now based on
the intersection of two instead of the union. Tests added.
This commit is contained in:
parent
934d091249
commit
eacc9b19a3
2 changed files with 36 additions and 9 deletions
|
@ -1472,8 +1472,8 @@ def _covariance_element(obs1, obs2):
|
|||
"""Estimates the covariance of two Obs objects, neglecting autocorrelations."""
|
||||
|
||||
def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx):
|
||||
deltas1 = _expand_deltas_for_merge(deltas1, idx1, len(idx1), new_idx)
|
||||
deltas2 = _expand_deltas_for_merge(deltas2, idx2, len(idx2), new_idx)
|
||||
deltas1 = _collapse_deltas_for_merge(deltas1, idx1, len(idx1), new_idx)
|
||||
deltas2 = _collapse_deltas_for_merge(deltas2, idx2, len(idx2), new_idx)
|
||||
return np.sum(deltas1 * deltas2)
|
||||
|
||||
if set(obs1.names).isdisjoint(set(obs2.names)):
|
||||
|
@ -1493,29 +1493,30 @@ def _covariance_element(obs1, obs2):
|
|||
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]])
|
||||
idl_d[r_name] = _intersection_idx([obs1.idl[r_name], obs2.idl[r_name]])
|
||||
|
||||
gamma = 0.0
|
||||
|
||||
for r_name in obs1.e_content[e_name]:
|
||||
if r_name not in obs2.e_content[e_name]:
|
||||
continue
|
||||
if len(idl_d[r_name]) == 0:
|
||||
continue
|
||||
gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name])
|
||||
|
||||
if gamma == 0.0:
|
||||
continue
|
||||
|
||||
gamma_div = 0.0
|
||||
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])
|
||||
e_N += len(idl_d[r_name])
|
||||
gamma /= max(gamma_div, 1.0)
|
||||
if len(idl_d[r_name]) == 0:
|
||||
continue
|
||||
gamma_div += np.sqrt(calc_gamma(obs1.deltas[r_name], obs1.deltas[r_name], obs1.idl[r_name], obs1.idl[r_name], idl_d[r_name]) * calc_gamma(obs2.deltas[r_name], obs2.deltas[r_name], obs2.idl[r_name], obs2.idl[r_name], idl_d[r_name]))
|
||||
gamma /= gamma_div
|
||||
|
||||
# Bias correction hep-lat/0306017 eq. (49)
|
||||
dvalue += (1 + 1 / e_N) * gamma / e_N
|
||||
dvalue += gamma
|
||||
|
||||
for e_name in obs1.cov_names:
|
||||
|
||||
|
|
|
@ -749,6 +749,32 @@ def test_covariance_idl():
|
|||
pe.covariance([obs1, obs2])
|
||||
|
||||
|
||||
def test_correlation_intersection_of_idls():
|
||||
range1 = range(1, 2000, 2)
|
||||
range2 = range(2, 2001, 2)
|
||||
|
||||
obs1 = pe.Obs([np.random.normal(1.0, 0.1, len(range1))], ["ens"], idl=[range1])
|
||||
obs2_a = 0.4 * pe.Obs([np.random.normal(1.0, 0.1, len(range1))], ["ens"], idl=[range1]) + 0.6 * obs1
|
||||
obs1.gamma_method()
|
||||
obs2_a.gamma_method()
|
||||
|
||||
cov1 = pe.covariance([obs1, obs2_a])
|
||||
corr1 = pe.covariance([obs1, obs2_a], correlation=True)
|
||||
|
||||
obs2_b = obs2_a + pe.Obs([np.random.normal(1.0, 0.1, len(range2))], ["ens"], idl=[range2])
|
||||
obs2_b.gamma_method()
|
||||
|
||||
cov2 = pe.covariance([obs1, obs2_b])
|
||||
corr2 = pe.covariance([obs1, obs2_b], correlation=True)
|
||||
|
||||
assert np.isclose(corr1[0, 1], corr2[0, 1], atol=1e-14)
|
||||
assert cov1[0, 1] > cov2[0, 1]
|
||||
|
||||
obs2_c = pe.Obs([np.random.normal(1.0, 0.1, len(range2))], ["ens"], idl=[range2])
|
||||
obs2_c.gamma_method()
|
||||
assert np.isclose(0, pe.covariance([obs1, obs2_c])[0, 1], atol=1e-14)
|
||||
|
||||
|
||||
def test_empty_obs():
|
||||
o = pe.Obs([np.random.rand(100)], ['test'])
|
||||
q = o + pe.Obs([], [], means=[])
|
||||
|
|
Loading…
Add table
Reference in a new issue