Adding variance of variance and normality tests for tally statistics (#3454)

Co-authored-by: Ethan Peterson <eepeterson3@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

tests/unit_tests/test_tallies.py

@@ -1,6 +1,8 @@
+from math import sqrt
 import numpy as np
 import pytest
 import openmc
+import scipy.stats as sps
 
 
 def test_xml_roundtrip(run_in_tmpdir):
@@ -163,3 +165,214 @@ def test_tally_application(sphere_model, run_in_tmpdir):
     assert (sp_tally.std_dev == tally.std_dev).all()
     assert (sp_tally.mean == tally.mean).all()
     assert sp_tally.nuclides == tally.nuclides
+
+def _tally_from_data(x, *, higher_moments=True, normality=True):
+    t = openmc.Tally()
+    t.scores = ["flux"]  # 1 score
+    t.nuclides = [openmc.Nuclide("H1")]  # 1 nuclide
+    t._sp_filename = "dummy.h5"  # mark "results available"
+    t._results_read = True  # don't try to read from disk
+    t._num_realizations = int(len(x))  # n
+    t.higher_moments = bool(higher_moments)
+
+    x = np.asarray(x, dtype=float)
+    # (num_filter_bins=1, num_nuclides=1, num_scores=1) -> (1,1,1) arrays
+    t._sum = np.array([[[np.sum(x)]]], dtype=float)
+    t._sum_sq = np.array([[[np.sum(x**2)]]], dtype=float)
+    if higher_moments:
+        t._sum_third = np.array([[[np.sum(x**3)]]], dtype=float)
+        t._sum_fourth = np.array([[[np.sum(x**4)]]], dtype=float)
+    return t
+
+@pytest.mark.parametrize(
+    "x, skew_true, kurt_true",
+    [   # Rademacher distribution
+        (np.array([1.0, -1.0] * 200), 0.0, 1.0),
+        # Two-point {0,3} with p(0)=3/4, p(3)=1/4
+        (np.concatenate([np.zeros(600), np.full(200, 3.0)]), 2.0 / sqrt(3.0), 7.0 / 3.0),
+        # Bernoulli distribution
+        (np.concatenate([np.ones(300), np.zeros(700)]), (1 - 2 * 0.3) / sqrt(0.3 * 0.7), (1 - 3 * 0.3 + 3 * 0.3**2) / (0.3 * 0.7)),
+    ],
+)
+def test_b1_b2_analytical_against_tally(x, skew_true, kurt_true):
+    t = _tally_from_data(x, higher_moments=True, normality=False)
+
+    g1 = t.skew(bias=True)[0, 0, 0]
+    b2 = t.kurtosis(bias=True, fisher=False)[0, 0, 0]
+
+    assert np.isclose(g1, skew_true, rtol=0, atol=1e-12)
+    assert np.isclose(b2, kurt_true, rtol=0, atol=1e-12)
+
+@pytest.mark.parametrize(
+    "draw, skew_true, kurt_true",
+    [(lambda rng, n: rng.normal(0, 1, n), 0.0, 3.0),  # Normal
+     (lambda rng, n: rng.random(n), 0.0, 1.8),  # Uniform(0,1)
+     (lambda rng, n: rng.exponential(1.0, n), 2.0, 9.0),  # Exp(1)
+     (lambda rng, n: (rng.random(n) < 0.3).astype(float),
+            (1 - 2 * 0.3) / sqrt(0.3 * 0.7),
+            (1 - 3 * 0.3 + 3 * 0.3**2) / (0.3 * 0.7),),],)
+
+def test_b1_b2_scipy_and_theory(draw, skew_true, kurt_true):
+    rng = np.random.default_rng(12345)
+    N = 200_000
+    x = draw(rng, N)
+
+    # Tally outputs
+    t = _tally_from_data(x, higher_moments=True, normality=False)
+    g1_t = t.skew(bias=True)[0, 0, 0]
+    b2_t = t.kurtosis(bias=True, fisher=False)[0, 0, 0]
+
+    # SciPy (population, bias=True to match population-moment style)
+    skew_sp = sps.skew(x, bias=True)
+    kurt_sp = sps.kurtosis(x, fisher=False, bias=True)
+
+    # Compare to SciPy numerically
+    assert np.isclose(g1_t, skew_sp, rtol=0, atol=5e-3)
+    assert np.isclose(b2_t, kurt_sp, rtol=0, atol=5e-3)
+
+    # Compare to analytical targets with size-dependent tolerances
+    tol_skew = 0.02 if abs(skew_true) < 0.5 else 0.05
+    tol_kurt = 0.03 if kurt_true < 4 else 0.1
+    assert abs(g1_t - skew_true) < tol_skew
+    assert abs(b2_t - kurt_true) < tol_kurt
+
+
+def test_kurtosis_bias_fisher_combinations():
+    """Test that all combinations of bias and fisher match scipy.stats.kurtosis"""
+    rng = np.random.default_rng(42)
+    x = rng.normal(0, 1, 10000)
+
+    t = _tally_from_data(x, higher_moments=True, normality=False)    # Test all four combinations
+    # 1. bias=True, fisher=False (Pearson's kurtosis, b2)
+    b2_tally = t.kurtosis(bias=True, fisher=False)[0, 0, 0]
+    b2_scipy = sps.kurtosis(x, fisher=False, bias=True)
+    assert np.isclose(b2_tally, b2_scipy, rtol=0, atol=1e-10)
+    assert np.isclose(b2_tally, 3.0, rtol=0.05, atol=0.1)  # Should be ~3 for normal
+
+    # 2. bias=True, fisher=True (excess kurtosis, g2)
+    g2_tally = t.kurtosis(bias=True, fisher=True)[0, 0, 0]
+    g2_scipy = sps.kurtosis(x, fisher=True, bias=True)
+    assert np.isclose(g2_tally, g2_scipy, rtol=0, atol=1e-10)
+    assert np.isclose(g2_tally, 0.0, rtol=0, atol=0.1)  # Should be ~0 for normal
+    assert np.isclose(g2_tally, b2_tally - 3.0, rtol=0, atol=1e-10)  # g2 = b2 - 3
+
+    # 3. bias=False, fisher=True (adjusted excess kurtosis, G2)
+    G2_tally = t.kurtosis(bias=False, fisher=True)[0, 0, 0]
+    G2_tally_default = t.kurtosis()[0, 0, 0]  # Should be same as default
+    G2_scipy = sps.kurtosis(x, fisher=True, bias=False)
+    assert np.isclose(G2_tally, G2_tally_default, rtol=0, atol=1e-10)
+    assert np.isclose(G2_tally, G2_scipy, rtol=0, atol=1e-10)
+    assert np.isclose(G2_tally, 0.0, rtol=0, atol=0.1)  # Should be ~0 for normal
+
+    # 4. bias=False, fisher=False (adjusted Pearson's kurtosis)
+    adj_b2_tally = t.kurtosis(bias=False, fisher=False)[0, 0, 0]
+    adj_b2_scipy = sps.kurtosis(x, fisher=False, bias=False)
+    assert np.isclose(adj_b2_tally, adj_b2_scipy, rtol=0, atol=1e-10)
+    assert np.isclose(adj_b2_tally, 3.0, rtol=0.05, atol=0.1)  # Should be ~3 for normal
+    assert np.isclose(adj_b2_tally, G2_tally + 3.0, rtol=0, atol=1e-10)  # adj_b2 = G2 + 3
+
+
+def test_ztests_scipy_comparison():
+    rng = np.random.default_rng(987)
+    x_norm = rng.normal(size=50_000)
+    x_exp = rng.exponential(size=50_000)
+
+    # -------- Normal dataset (should not reject) --------
+    t0 = _tally_from_data(x_norm, higher_moments=True, normality=True)
+    Zb1_0, p_skew_0 = t0.skewtest(alternative="two-sided")
+    Zb2_0, p_kurt_0 = t0.kurtosistest(alternative="two-sided")
+    K2_0, p_omni_0 = t0.normaltest(alternative="two-sided")
+
+    Zb1_0 = Zb1_0.ravel()[0]
+    p_skew_0 = p_skew_0.ravel()[0]
+    Zb2_0 = Zb2_0.ravel()[0]
+    p_kurt_0 = p_kurt_0.ravel()[0]
+    K2_0 = K2_0.ravel()[0]
+    p_omni_0 = p_omni_0.ravel()[0]
+
+    z_skew_sp0, p_skew_sp0 = sps.skewtest(x_norm)
+    z_kurt_sp0, p_kurt_sp0 = sps.kurtosistest(x_norm)
+    k2_sp0, p_omni_sp0 = sps.normaltest(x_norm)
+
+    assert np.isclose(Zb1_0, z_skew_sp0, atol=0.15)
+    assert np.isclose(Zb2_0, z_kurt_sp0, atol=0.15)
+    assert np.isclose(K2_0, k2_sp0, atol=0.30)
+    assert np.isclose(p_skew_0, p_skew_sp0, atol=5e-3)
+    assert np.isclose(p_kurt_0, p_kurt_sp0, atol=5e-3)
+    assert np.isclose(p_omni_0, p_omni_sp0, atol=5e-3)
+
+    # -------- Exponential dataset (should strongly reject) --------
+    t1 = _tally_from_data(x_exp, higher_moments=True, normality=True)
+
+    Zb1_1, p_skew_1 = t1.skewtest(alternative="two-sided")
+    Zb2_1, p_kurt_1 = t1.kurtosistest(alternative="two-sided")
+    K2_1, p_omni_1 = t1.normaltest(alternative="two-sided")
+
+    Zb1_1 = Zb1_1.ravel()[0]
+    p_skew_1 = p_skew_1.ravel()[0]
+    Zb2_1 = Zb2_1.ravel()[0]
+    p_kurt_1 = p_kurt_1.ravel()[0]
+    K2_1 = K2_1.ravel()[0]
+    p_omni_1 = p_omni_1.ravel()[0]
+
+    z_skew_sp1, p_skew_sp1 = sps.skewtest(x_exp)
+    z_kurt_sp1, p_kurt_sp1 = sps.kurtosistest(x_exp)
+    k2_sp1, p_omni_sp1 = sps.normaltest(x_exp)
+
+    # Both pipelines should reject very strongly
+    assert p_skew_1 < 1e-6 and p_skew_sp1 < 1e-6
+    assert p_kurt_1 < 1e-6 and p_kurt_sp1 < 1e-6
+    assert p_omni_1 < 1e-6 and p_omni_sp1 < 1e-6
+
+    # Right-skewed and heavy-tailed → large positive Z-statistics
+    assert Zb1_1 > 30 and z_skew_sp1 > 30
+    assert Zb2_1 > 30 and z_kurt_sp1 > 30
+    assert K2_1 > 2000 and k2_sp1 > 2000
+
+def test_vov_stochastic(sphere_model, run_in_tmpdir):
+    tally = openmc.Tally(name="test tally")
+    ef = openmc.EnergyFilter([0.0, 0.1, 1.0, 10.0e6])
+    mesh = openmc.RegularMesh.from_domain(sphere_model.geometry, (2, 2, 2))
+    mf = openmc.MeshFilter(mesh)
+    tally.filters = [ef, mf]
+    tally.scores = ["flux", "absorption", "fission", "scatter"]
+    tally.higher_moments = True
+    sphere_model.tallies = [tally]
+
+    sp_file = sphere_model.run(apply_tally_results=True)
+
+    assert tally._mean is None
+    assert tally._std_dev is None
+    assert tally._sum is None
+    assert tally._sum_sq is None
+    assert tally._sum_third is None
+    assert tally._sum_fourth is None
+    assert tally._num_realizations == 0
+    assert tally._sp_filename == sp_file
+
+    with openmc.StatePoint(sp_file) as sp:
+        assert tally in sp.tallies.values()
+        sp_tally = sp.tallies[tally.id]
+
+    assert np.all(sp_tally.std_dev == tally.std_dev)
+    assert np.all(sp_tally.mean == tally.mean)
+    assert np.all(sp_tally.vov == tally.vov)
+    assert sp_tally.nuclides == tally.nuclides
+
+    n = sp_tally.num_realizations
+    mean = sp_tally.mean
+    sum_ = sp_tally._sum
+    sum_sq = sp_tally._sum_sq
+    sum_third = sp_tally._sum_third
+    sum_fourth = sp_tally._sum_fourth
+
+    expected_vov = np.zeros_like(mean)
+    nonzero = np.abs(mean) > 0
+
+    num = (sum_fourth - (4.0*sum_third*sum_)/n + (6.0*sum_sq*sum_**2)/(n**2)
+           - (3.0*sum_**4)/(n**3))
+    den = (sum_sq - (1.0/n)*sum_**2)**2
+
+    expected_vov[nonzero] = num[nonzero]/den[nonzero] - 1.0/n
+
+    assert np.allclose(expected_vov, sp_tally.vov, rtol=1e-7, atol=0.0)