Adding variance of variance and normality tests for tally statistics

docs/source/methods/tallies.rst

@@ -387,6 +387,98 @@ 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`) and 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`.
+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 [Agostino]_. 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 ``vov_enabled = True`` for 
+   running these diagnostics including the skewness, kurtosis and normality tests.
+
 Figure of Merit
 +++++++++++++++
 
@@ -412,7 +504,9 @@ 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 +615,9 @@ improve the estimate of the percentile.
 
    .. rubric:: References
 
+.. [Agostino] D'Agostino, Ralph B., and Albert Belanger. *A Suggestion for Using Powerful and Informative Tests of Normality.* 
+              The American Statistician, vol. 44, no. 4, (1990).
+
 .. _following approximation: https://doi.org/10.1080/03610918708812641
 
 .. _Bessel's correction: https://en.wikipedia.org/wiki/Bessel's_correction
@@ -541,6 +638,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

@@ -286,7 +286,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 vov_enabled() const { return vov_enabled_; }
+
   //----------------------------------------------------------------------------
   // Other methods.
 
@@ -190,6 +192,8 @@ class Tally {
   //! Whether to multiply by atom density for reaction rates
   bool multiply_density_ {true};
 
+  bool vov_enabled_ = {false};
+
   int64_t index_;
 };
 

openmc/statepoint.py

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