Adding variance of variance and normality tests for tally statistics

[Grego01-biot]

Description

This PR consists of extending the statistical capabilities of tallies in OpenMC by adding the variance of the variance, the skewness and kurtosis values (with adjusted moments) and the D'Agostino-Pearson normality tests (skewness, kurtosis, and omnibus) in fixed source and eigenvalue simulations [1],[2].

The variance of the variance involves the estimated third and fourth moments of the output distribution and is much more sensitive to large fluctuations. It is used to measure the relative statistical uncertainty in the estimated relative error. Since the third and fourth moments are accumulated, additional information can be extracted and normality tests can be performed [3].

The normality tests allow to quantify the nature of the non-normality to see if the data distribution is skewed to the right/left (skewness) and heavy/light in the tails (kurtosis). The most commonly used normality tests are tests against the null hypothesis H0 that the data is normally distributed. The first step in the testing procedure is to compute a test statistic, which follows some known distribution. The test statistic is then converted into a p-value. The null hypothesis is rejected if the p-value is less than a predetermined significance level that can be defined by the user in the tally settings. If the p-value is greater than the significance level, one fails to reject the null hypothesis. For example, if the significance level is chosen to be 0.05 and the obtained p-value is less than this, one can conclude that with 95% probability the data isn’t normally distributed [4].

Future work should focus on adding the capability to print how a tally mean, standard deviation, figure of merit and variance of the variance fluctuate as a function of the number of batches similar to the TFC in MCNP.

[1] D. N. Joanes and C. A. Gill. "Comparing Measures of Sample Skewness and Kurtosis". Journal of the Royal Statistical Society. Series D (The Statistician), 1998, Vol. 47, No. 1 (1998), pp. 183-189.
[2] 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, pp. 316–21. JSTOR, https://doi.org/10.2307/2684359. Accessed 17 June 2025.
[3] Kulesza, Joel A. and Adams, Terry R. and Armstrong, Jerawan C. and Bolding, Simon R. and Brown, Forrest B. and Bull, Jeffrey S. and Burke, Timothy P. and Clark, Alexander R. and Forster, III, Robert Arthur and Giron, Jesse F. and Grieve, Avery S. and Josey, Colin J. and Martz, Roger L. and McKinney, Gregg W. and Pearson, Eric J. and Rising, Michael E. and Solomon, Jr., Clell J. and Swaminarayan, Sriram and Trahan, Travis J. and Weaver, Colin A. and Wilson, Stephen C. and Zukaitis, Anthony J. "MCNP Code Version 6.3.1 Theory and User Manual", LA-UR-24-24602, Rev.~1, LANL, NM, USA, May 2024.
[4] T. Kaltiaisenaho. "Statistical Tests and the Underestimation of Variance in Serpent 2", Tech Report: VTT-R-00371-14: https://serpent.vtt.fi/serpent/download/VTT-R-00371-14.pdf

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[paulromano]

@Grego01-biot I went through and did a bit of refactoring here:

• Made m2, m3, and m4 properties
• Made skew and kurtosis methods that follow the same conventions as in scipy.stats
• Renamed vov_enabled to higher_moments

Please see my comments below; once addressed I think this will be very close to being merge-ready!

[paulromano]

Just finished going through this branch and making some changes. The changes I've made are mostly to have an API that is consistent with scipy when it comes to computing the skewness, kurtosis, and performing statistical tests. I've also removed the print_normality_test for now as we don't typically have methods in OpenMC that print lots of information like that -- I think that is best left as user code. @Grego01-biot If you want, you can put together a Jupyter notebook showing the user of the higher moments feature and we can add it to our list of example notebooks.

[eepeterson]

thanks @Grego01-biot and @paulromano for this highly asked for feature!