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>

Author: 3 people

docs/source/io_formats/statepoint.rst

@@ -149,6 +149,8 @@ The current version of the statepoint file format is 18.1.
                tallies will have a value of 0 unless otherwise instructed.
              - **multiply_density** (*int*) -- Flag indicating whether reaction
                rates should be multiplied by atom density (1) or not (0).
+             - **higher_moments** (*int*) -- Flag indicating whether
+               higher-order tally moments are enabled (1) or not (0).
 
 :Datasets: - **n_realizations** (*int*) -- Number of realizations.
            - **n_filters** (*int*) -- Number of filters used.

docs/source/methods/tallies.rst

@@ -387,6 +387,101 @@ of this is that the longer you run a simulation, the better you know your
 results. Therefore, by running a simulation long enough, it is possible to
 reduce the stochastic uncertainty to arbitrarily low levels.
 
+Skewness
+++++++++
+
+The `skewness`_ of a population quantifies the asymmetry of the probability
+distribution around its mean. Positive and negative skewness indicate a
+longer/heavier right and left tail respectively. Let :math:`x_1,\ldots,x_n` be
+the per-realization values for a bin, with sample mean :math:`\bar{x}` and
+sample central moments:
+
+.. math::
+
+   m_k \;=\; \frac{1}{n}\sum_{i=1}^{n}\bigl(x_i-\bar{x}\bigr)^k.
+
+OpenMC reports the *adjusted Fisher-Pearson skewness* (defined for :math:`n \ge
+3`), which is commonly used in many statistical packages:
+
+.. math::
+
+   G_1 \;=\; \frac{\sqrt{n \cdot (n-1)}}{\,n-2\,}\cdot\frac{m_3}{m_2^{3/2}}.
+
+where :math:`m_2` and :math:`m_3` correspond to the biased sample second and
+third central moment respectively.
+
+Kurtosis
+++++++++
+
+The `kurtosis`_ of a population quantifies tail weight (also called tailedness)
+of the probability distribution relative to a normal distribution. Positive
+excess kurtosis indicates *heavier tails* whereas negative excess kurtosis
+indicates *lighter tails*. Kurtosis is especially useful for identifying bins
+where occasional extreme scores dominate uncertainty. OpenMC reports the
+*adjusted excess kurtosis* (defined for :math:`n \ge 4`):
+
+.. math::
+
+   G_2 \;=\; \frac{(n-1)}{(n-2)(n-3)}
+             \left[(n+1)\,\frac{m_4}{m_2^{2}} \;-\; 3(n-1)\right].
+
+where :math:`m_2` and :math:`m_4` correspond to the biased sample second and
+fourth central moment respectively. For a perfectly normal distribution, the
+excess kurtosis is :math:`0`.
+
+Variance of Variance
+++++++++++++++++++++
+
+The variance of the variance (also known as the coefficient of variation
+squared) measures *stability of the sample variance* :math:`s^2` and, by
+extension, the reliability of reported relative errors. High VOV means that
+error bars themselves are noisy—often due to heavy tails, skewness, or too few
+realizations.
+
+.. math::
+
+    VOV = \frac{s^2(s_{\bar{X}}^2)}{s_{\bar{X}}^4 } = \frac{m_4}{m_2^2} - \frac{1}{n}
+
+where :math:`s_{\bar{X}}^2` is the estimated variance of the mean and
+:math:`s^2(s_{\bar{X}}^2)` is the estimated variance in :math:`s_{\bar{X}}^2`.
+The MCNP manual suggests a hard threshold such that :math:`VOV < 0.1` to improve
+the probability of forming a reliable confidence interval. However, OpenMC does
+not enforce an universal cut-off because the suitability of any single threshold
+depends strongly on problem specifics (estimator choice, variance-reduction
+settings, tally binning, or even effective sample size).
+
+
+Normality Tests (D'Agostino-Pearson)
+++++++++++++++++++++++++++++++++++++
+
+These normality test verify the hypothesis that fluctuations are *approximately
+normal*, a working assumption behind many Monte Carlo diagnostics and
+`confidence-interval heuristics`_. Tests are provided for: (i) skewness-only,
+(ii) kurtosis-only, and (iii) the *omnibus* combination. OpenMC uses the
+finite-sample-adjusted skewness :math:`G_1` and excess kurtosis :math:`G_2`
+above to construct standardized normal scores :math:`Z_1` (from :math:`G_1`) and
+:math:`Z_2` (from :math:`G_2`) via the D'Agostino-Pearson transformations. The
+omnibus statistic is
+
+.. math::
+
+   K^2 \;=\; Z_1^{\,2} \;+\; Z_2^{\,2}
+   \;\sim\; \chi^2_{(2)} \quad \text{under } H_0:\ \text{normality}.
+
+OpenMC reports :math:`Z_1`, :math:`Z_2`, :math:`K^2`, and their p-values when
+prerequisites are met (skewness for :math:`n\ge 3`, kurtosis and omnibus for
+:math:`n\ge 4`). Given a user-chosen significance level :math:`\alpha` (default
+is :math:`0.05`), reject :math:`H_0` if :math:`\text{p-value}<\alpha`; otherwise
+fail to reject. OpenMC leaves the interpretation to the user, who should
+consider VOV together with skewness, kurtosis, and normality tests results when
+judging whether reported confidence intervals are credible for their application
+[#norm-tests]_.
+
+.. [#norm-tests]
+   Higher-moments accumulation must be enabled with ``higher_moments = True``
+   for running these diagnostics including the skewness, kurtosis, and normality
+   tests.
+
 Figure of Merit
 +++++++++++++++
 
@@ -405,14 +500,16 @@ defined as
 .. math::
     :label: relative_error
 
-    r = \frac{s_\bar{X}}{\bar{x}}.
+    r = \frac{s_{\bar{X}}}{\bar{x}}.
 
 Based on this definition, one can see that a higher FOM is desirable. The FOM is
 useful as a comparative tool. For example, if a variance reduction technique is
 being applied to a simulation, the FOM with variance reduction can be compared
 to the FOM without variance reduction to ascertain whether the reduction in
 variance outweighs the potential increase in execution time (e.g., due to
-particle splitting).
+particle splitting). It is important to note that MCNP reports the FOM using CPU
+time (wall-clock time multiplied by the number of threads/cores), whereas OpenMC
+reports the FOM using only the wall-clock time :math:`t`.
 
 Confidence Intervals
 ++++++++++++++++++++
@@ -521,6 +618,8 @@ improve the estimate of the percentile.
 
    .. rubric:: References
 
+.. _confidence-interval heuristics: https://doi.org/10.1080/00031305.1990.10475751
+
 .. _following approximation: https://doi.org/10.1080/03610918708812641
 
 .. _Bessel's correction: https://en.wikipedia.org/wiki/Bessel's_correction
@@ -541,6 +640,10 @@ improve the estimate of the percentile.
 
 .. _converges in distribution: https://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_distribution
 
+.. _skewness: https://en.wikipedia.org/wiki/Skewness
+
+.. _kurtosis: https://en.wikipedia.org/wiki/Kurtosis
+
 .. _confidence intervals: https://en.wikipedia.org/wiki/Confidence_interval
 
 .. _Student's t-distribution: https://en.wikipedia.org/wiki/Student%27s_t-distribution

include/openmc/constants.h

@@ -291,7 +291,7 @@ enum class MgxsType {
 // ============================================================================
 // TALLY-RELATED CONSTANTS
 
-enum class TallyResult { VALUE, SUM, SUM_SQ, SIZE };
+enum class TallyResult { VALUE, SUM, SUM_SQ, SUM_THIRD, SUM_FOURTH };
 
 enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };
 

include/openmc/hdf5_interface.h

@@ -100,8 +100,8 @@ void read_llong(hid_t obj_id, const char* name, long long* buffer, bool indep);
 void read_string(
   hid_t obj_id, const char* name, size_t slen, char* buffer, bool indep);
 
-void read_tally_results(
-  hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results);
+void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
+  hsize_t n_results, double* results);
 void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims,
   const char* name, const double* buffer);
 void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims,
@@ -114,9 +114,9 @@ void write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
 void write_llong(hid_t group_id, int ndim, const hsize_t* dims,
   const char* name, const long long* buffer, bool indep);
 void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
-  const char* name, char const* buffer, bool indep);
-void write_tally_results(
-  hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results);
+  const char* name, const char* buffer, bool indep);
+void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
+  hsize_t n_results, const double* results);
 } // extern "C"
 
 //==============================================================================

include/openmc/tallies/tally.h

@@ -106,6 +106,8 @@ class Tally {
 
   bool writable() const { return writable_; }
 
+  bool higher_moments() const { return higher_moments_; }
+
   //----------------------------------------------------------------------------
   // Other methods.
 
@@ -190,6 +192,9 @@ class Tally {
   //! Whether to multiply by atom density for reaction rates
   bool multiply_density_ {true};
 
+  //! Whether to accumulate higher moments (third and fourth)
+  bool higher_moments_ {false};
+
   int64_t index_;
 };
 

openmc/statepoint.py

@@ -434,6 +434,10 @@ def tallies(self):
                     if "multiply_density" in group.attrs:
                         tally.multiply_density = group.attrs["multiply_density"].item() > 0
 
+                    # Check if tally has higher_moments attribute
+                    if 'higher_moments' in group.attrs:
+                        tally.higher_moments = bool(group.attrs['higher_moments'][()])
+
                     # Read the number of realizations
                     n_realizations = group['n_realizations'][()]