Hypothesis testing of scientific Monte Carlo calculations

The steadily increasing size of scientific Monte Carlo simulations and the desire for robust, correct, and reproducible results necessitates rigorous testing procedures for scientific simulations in order to detect numerical problems and programming bugs. However, the testing paradigms developed for deterministic algorithms have proven to be ill suited for stochastic algorithms. In this paper we demonstrate explicitly how the technique of statistical hypothesis testing, which is in wide use in other fields of science, can be used to devise automatic and reliable tests for Monte Carlo methods, and we show that these tests are able to detect some of the common problems encountered in stochastic scientific simulations. We argue that hypothesis testing should become part of the standard testing toolkit for scientific simulations.

Figure 1

Scalar one-sample test against $\langle m\rangle =0$ for single spin-flip updates (red curves) and Wolff cluster updates (blue curves) in a classical two-dimensional Ising model with length $L=16$: (a) result for $\langle m\rangle$ from $N^{\prime }={10}^6$ Monte Carlo sweeps (a sweep is either a set of $L^2$ single spin flips or a single cluster update); (b) $\left|t\right|$ score as the left-hand side of Eq. (2); and (c) $p$ values from a two-tailed test of the $t$ score against the Student $t_{N-1}$ distribution (the shaded area indicates $p<\alpha =0.01$ and thus a failed test).

Figure 2

Scalar two-sample test for the Binder cumulant $U_4$ in the Ising model obtained from Markov chain Monte Carlo with Wolff cluster updates at temperature $T=2.3$ for system length $L=32$. (a) Sample mean of $U_4$ of an erroneous implementation for different simulation times $N$, and the result of a correct reference simulation; (b) corresponding $p$ values computed using Eq. (6); (c) illustration of the $p$ value for $N^{\prime }={10}^3$ $(N=123$ after the removal of autocorrelation and variance pooling) as the shaded area under the probability density function (p.d.f.) of the corresponding $t$ distribution.

Figure 3

Vector-valued two-sample test on the spin correlation function ${\chi }_{x,0}$ [cf. Eq. (11)]. (a) Simulation result for Wolff updates (black curve) and single spin-flip updates (red curve); (b) deviation of spin-flip from Wolff update (the shaded region are the Wolff result error bars); (c) projected deviations and errors (numerator and denominator in the l.h.s. of Eq. (14)); (d) correlation matrix (12); (e) p.d.f. of the corresponding $F$ distribution in Eq. (14) with the mean as vertical blue line and $p_{<}$ value as blue shaded area to the left of the observed $T^2$ score (blue dot) as well as test failure threshold as red-shaded area.

Figure 4

(a) Green's function $G\left(\tau \right)$ for AIM parameters as in main text: Monte Carlo result (red) and exact result with an artificially introduced shift of the half bin size, i.e., $G(\tau -0.02)$, (black) modeling a binning error; (b) deviation from the exact result in multiples of the standard error as well as a linear regression (black dashed line) $0.539-0.098\tau$ with a goodness of fit of $R^2\approx 0.065$.