A tail-regression estimator for heavy-tailed distributions of known tail indices and its application to continuum quantum Monte Carlo data

TL;DR

This paper introduces a regression-based estimator for heavy-tailed distributions with known tail indices, enabling reliable confidence intervals in quantum Monte Carlo data analysis where traditional methods fail.

Contribution

It presents a novel tail-regression estimator that leverages known asymptotic decay exponents to accurately estimate moments and confidence intervals for heavy-tailed distributions.

Findings

Accurately estimates variance of local energy in electron gas

Provides confidence intervals for Hellmann-Feynman forces in carbon dimer

Reduces uncertainty by 45-60 times compared to standard error

Abstract

Standard statistical analysis is unable to provide reliable confidence intervals on expectation values of probability distributions that do not satisfy the conditions of the central limit theorem. We present a regression-based estimator of an arbitrary moment of a probability distribution with power-law heavy tails that exploits knowledge of the exponents of its asymptotic decay to bypass this issue entirely. Our method is applied to synthetic data and to energy and atomic force data from variational and diffusion quantum Monte Carlo calculations, whose distributions have known asymptotic forms [J. R. Trail, Phys. Rev. E 77, 016703 (2008); A. Badinski et al., J. Phys.: Condens. Matter 22 074202 (2010)]. We obtain convergent, accurate confidence intervals on the variance of the local energy of an electron gas and on the Hellmann-Feynman force on an atom in the all-electron carbon dimer. In each of these cases the uncertainty on our estimator is 45% and 60 times smaller, respectively, than the nominal (ill-defined) standard error.

Tables

Table 1. Table 1: Standard and tail-regression estimators of ⟨ A ⟩ delimited-⟨⟩ 𝐴 \langle A\rangle and σ 2 superscript 𝜎 2 \sigma^{2} for model distribution P A ​ ( A ) = 1 2 ​ H 3.1 ​ ( A ) + 1 2 ​ H 4.1 ​ ( A ) subscript 𝑃 𝐴 𝐴 1 2 subscript 𝐻 3.1 𝐴 1 2 subscript 𝐻 4.1 𝐴 P_{A}(A)=\frac{1}{2}H_{3.1}(A)+\frac{1}{2}H_{4.1}(A) obtained from random samples of various sizes M 𝑀 M . The “optimal” expansion orders n R subscript 𝑛 R n_{\rm R} and thresholds − log ⁡ q R subscript 𝑞 R -\log q_{\rm R} used for the tail-regression estimator in each case are also shown.

		Standard					TRE (unconstrained)							TRE (1 constraint)
$M$		$\overline{A}$		$S^2$			$n_R$	$-\log q_R$	$𝒜$		$𝒱$			$n_R$	$-\log q_R$	$𝒜$		$𝒱$
${10}^3$		$0$.	$001\left(39\right)$	$1$.	$55\left(28\right)$		$1$	$3.25$	$0$.	$006\left(38\right)$	$3$.	$8\left(17\right)$		$5$	$1.00$	$0$.	$001\left(34\right)$	$3$.	$6\left(20\right)$
${10}^4$		$-0$.	$000\left(18\right)$	$3$.	$3\left(13\right)$		$4$	$1.25$	$-0$.	$000\left(14\right)$	$2$.	$71\left(73\right)$		$3$	$1.50$	$0$.	$001\left(11\right)$	$3$.	$18\left(65\right)$
${10}^5$		$0$.	$0009\left(50\right)$	$2$.	$51\left(31\right)$		$4$	$1.50$	$0$.	$0009\left(48\right)$	$4$.	$54\left(33\right)$		$4$	$1.50$	$0$.	$0057\left(35\right)$	$4$.	$55\left(33\right)$
${10}^6$		$-0$.	$0017\left(16\right)$	$2$.	$63\left(19\right)$		$3$	$2.25$	$-0$.	$0014\left(15\right)$	$4$.	$57\left(12\right)$		$3$	$2.25$	$-0$.	$0011\left(11\right)$	$4$.	$56\left(12\right)$
${10}^7$		$-0$.	$00022\left(52\right)$	$2$.	$69\left(12\right)$		$2$	$3.75$	$-0$.	$00017\left(46\right)$	$4$.	$585\left(41\right)$		$3$	$2.75$	$-0$.	$00032\left(37\right)$	$4$.	$608\left(49\right)$
${10}^8$		$0$.	$00018\left(16\right)$	$2$.	$689\left(66\right)$		$3$	$2.75$	$0$.	$00016\left(15\right)$	$4$.	$661\left(15\right)$		$3$	$2.75$	$0$.	$00019\left(12\right)$	$4$.	$662\left(15\right)$
Exact		$0$.	$0$	$4$.	$6586$				$0$.	$0$	$4$.	$6586$				$0$.	$0$	$4$.	$6586$

Table 2. Table 2: Standard and tail-regression estimators of ⟨ A ⟩ delimited-⟨⟩ 𝐴 \langle A\rangle for model distribution P A ​ ( A ) = 1 2 ​ H 2.1 ​ ( A ) + 1 2 ​ H 3.1 ​ ( A ) subscript 𝑃 𝐴 𝐴 1 2 subscript 𝐻 2.1 𝐴 1 2 subscript 𝐻 3.1 𝐴 P_{A}(A)=\frac{1}{2}H_{2.1}(A)+\frac{1}{2}H_{3.1}(A) obtained from random samples of various sizes M 𝑀 M . The “optimal” expansion orders n R subscript 𝑛 R n_{\rm R} and thresholds − log ⁡ q R subscript 𝑞 R -\log q_{\rm R} used for the tail-regression estimator in each case are also shown.

		Standard			TRE (unconstrained)					TRE (1 constraint)
$M$		$\overline{A}$			$n_R$	$-\log q_R$	$𝒜$			$n_R$	$-\log q_R$	$𝒜$
${10}^3$		$-6$.	$4\left(67\right)$		$1$	$1.30$	$-0$.	$12\left(29\right)$		$1$	$1.60$	$-0$.	$026\left(43\right)$
${10}^4$		$-0$.	$69\left(69\right)$		$5$	$1.00$	$-0$.	$11\left(20\right)$		$4$	$1.00$	$-0$.	$005\left(18\right)$
${10}^5$		$0$.	$18\left(19\right)$		$5$	$1.10$	$0$.	$005\left(72\right)$		$6$	$1.00$	$-0$.	$0131\left(76\right)$
${10}^6$		$-0$.	$21\left(17\right)$		$7$	$1.00$	$-0$.	$007\left(26\right)$		$4$	$1.40$	$-0$.	$0021\left(25\right)$
${10}^7$		$-0$.	$115\left(66\right)$		$5$	$1.30$	$0$.	$0084\left(84\right)$		$4$	$1.50$	$0$.	$00073\left(83\right)$
${10}^8$		$-0$.	$102\left(66\right)$		$4$	$1.70$	$0$.	$0028\left(28\right)$		$4$	$1.70$	$-0$.	$00003\left(29\right)$
Exact		$0$.	$0$				$0$.	$0$				$0$.	$0$

Table 3. Table 3: Standard and tail-regression estimators of ⟨ A ⟩ delimited-⟨⟩ 𝐴 \langle A\rangle for model distribution P A ​ ( A ) = 1 2 ​ H 1.1 ​ ( A ) + 1 2 ​ H 2.1 ​ ( A ) subscript 𝑃 𝐴 𝐴 1 2 subscript 𝐻 1.1 𝐴 1 2 subscript 𝐻 2.1 𝐴 P_{A}(A)=\frac{1}{2}H_{1.1}(A)+\frac{1}{2}H_{2.1}(A) obtained from random samples of various sizes M 𝑀 M . The “optimal” expansion orders n R subscript 𝑛 R n_{\rm R} and thresholds − log ⁡ q R subscript 𝑞 R -\log q_{\rm R} used for the tail-regression estimator in each case are also shown.

		Standard		TRE (1 constraint)
$M$		$\left|\overline{A}\right|$		$n_R$	$-\log q_R$	$𝒜$
${10}^3$		$\sim {10}^{30}$		$1$	$1.15$	$-0$.	$42\left(55\right)$
${10}^4$		$\sim {10}^{28}$		$2$	$0.95$	$0$.	$01\left(18\right)$
${10}^5$		$\sim {10}^{39}$		$2$	$1.00$	$0$.	$046\left(68\right)$
${10}^6$		$\sim {10}^{48}$		$3$	$1.10$	$-0$.	$028\left(54\right)$
${10}^7$		$\sim {10}^{59}$		$6$	$0.85$	$-0$.	$019\left(20\right)$
${10}^8$		$\sim {10}^{79}$		$5$	$0.95$	$-0$.	$0062\left(76\right)$
Exact		$0.0$				$0$.	$0$

Table 4. Table 4: Standard and tail-regression estimators of ⟨ A ⟩ delimited-⟨⟩ 𝐴 \langle A\rangle and σ 2 superscript 𝜎 2 \sigma^{2} for the VMC energy of the 54-electron gas at r s = 1 subscript 𝑟 s 1 r_{\rm s}=1 using the SJ wave function obtained from local energy samples of various sizes M 𝑀 M . The “optimal” expansion order n R subscript 𝑛 R n_{\rm R} and threshold − log ⁡ q R subscript 𝑞 R -\log q_{\rm R} used for the tail-regression estimator in each case are also shown.

		Standard					TRE (unconstrained)							TRE (1 constraint)
$M$		$\overline{A}$		$S^2$			$n_R$	$-\log q_R$	$𝒜$		$𝒱$			$n_R$	$-\log q_R$	$𝒜$		$𝒱$
${10}^3$		$0$.	$53381\left(74\right)$	$0$.	$000542\left(34\right)$		$2$	$4.25$	$0$.	$5339\left(12\right)$	$0$.	$00090\left(30\right)$		$2$	$4.50$	$0$.	$53342\left(96\right)$	$0$.	$00095\left(17\right)$
${10}^4$		$0$.	$53314\left(24\right)$	$0$.	$000581\left(13\right)$		$2$	$6.50$	$0$.	$53309\left(26\right)$	$0$.	$000646\left(40\right)$		$1$	$7.75$	$0$.	$53310\left(26\right)$	$0$.	$000627\left(28\right)$
${10}^5$		$0$.	$532868\left(76\right)$	$0$.	$0005738\left(39\right)$		$1$	$8.00$	$0$.	$532863\left(78\right)$	$0$.	$0005843\left(61\right)$		$1$	$8.25$	$0$.	$532868\left(76\right)$	$0$.	$0005839\left(57\right)$
${10}^6$		$0$.	$532820\left(24\right)$	$0$.	$0005804\left(44\right)$		$2$	$7.25$	$0$.	$532817\left(24\right)$	$0$.	$0005806\left(24\right)$		$2$	$7.25$	$0$.	$532815\left(24\right)$	$0$.	$0005805\left(24\right)$
${10}^7$		$0$.	$5328407\left(76\right)$	$0$.	$0005761\left(15\right)$		$2$	$8.00$	$0$.	$5328396\left(76\right)$	$0$.	$00057618\left(82\right)$		$3$	$8.25$	$0$.	$5328396\left(76\right)$	$0$.	$0005763\left(12\right)$
${10}^8$		$0$.	$5328500\left(24\right)$	$0$.	$00057442\left(31\right)$		$2$	$8.25$	$0$.	$5328500\left(24\right)$	$0$.	$00057545\left(25\right)$		$2$	$8.25$	$0$.	$5328499\left(24\right)$	$0$.	$00057545\left(24\right)$

Table 5. Table 5: Standard and tail-regression estimators of ⟨ A ⟩ delimited-⟨⟩ 𝐴 \langle A\rangle for the VMC Hellmann-Feynman force on a carbon atom along the molecular axis of the C 2 molecule at r CC = 2 subscript 𝑟 CC 2 r_{\rm CC}=2 a.u., obtained from local force samples of various sizes M 𝑀 M . The “optimal” expansion orders n R subscript 𝑛 R n_{\rm R} and thresholds − log ⁡ q R subscript 𝑞 R -\log q_{\rm R} used for the tail-regression estimator in each case are also shown. The “best” value is provided for reference and corresponds to the tail-regression estimator of ⟨ F HFT + F ZV ⟩ delimited-⟨⟩ subscript 𝐹 HFT subscript 𝐹 ZV \langle F_{\rm HFT}+F_{\rm ZV}\rangle using 10 8 superscript 10 8 10^{8} sample points.

		Standard			TRE (unconstrained)					TRE (1 constraint)
$M$		$\overline{A}$			$n_R$	$-\log q_R$	$𝒜$			$n_R$	$-\log q_R$	$𝒜$
${10}^3$		$-26$	$\left(56\right)$		$4$	$0.90$	$-18$	$\left(24\right)$		$3$	$0.90$	$0$	$.7\left(78\right)$
${10}^4$		$-29$	$\left(22\right)$		$5$	$0.80$	$-3$	$.3\left(88\right)$		$4$	$0.90$	$-0$	$.7\left(33\right)$
${10}^5$		$-2$	$.1\left(80\right)$		$5$	$0.80$	$0$	$.3\left(28\right)$		$5$	$0.80$	$0$	$.3\left(10\right)$
${10}^6$		$-17$	$.0\left(73\right)$		$4$	$1.00$	$-0$	$.1\left(10\right)$		$4$	$1.00$	$-0$	$.34\left(40\right)$
${10}^7$		$-8$	$.0\left(79\right)$		$3$	$1.30$	$-0$	$.43\left(34\right)$		$3$	$1.30$	$-0$	$.50\left(13\right)$
${10}^8$		$-0$	$.3\left(15\right)$		$3$	$1.40$	$-0$	$.51\left(11\right)$		$3$	$1.40$	$-0$	$.567\left(43\right)$
Best		$-0$	$.5770\left(29\right)$				$-0$	$.5770\left(29\right)$				$-0$	$.5770\left(29\right)$

Table 6. Table 6: Standard and tail-regression estimators of ⟨ A ⟩ delimited-⟨⟩ 𝐴 \langle A\rangle for the VMC local Hellmann-Feynman force, zero-variance corrected Hellmann-Feynman force, Pulay force, and zero-variance corrected total force on a carbon atom along the molecular axis of the C 2 molecule at r CC = 2 subscript 𝑟 CC 2 r_{\rm CC}=2 a.u., obtained from 10 8 superscript 10 8 10^{8} sample points.

	Standard		TRE		TRE(1)
$⟨F_{\mathrm{HFT}}⟩$	$-0$	$.3\left(15\right)$	$-0$	$.51\left(11\right)$	$-0$	$.567\left(43\right)$
$⟨F_{\mathrm{HFT}}+F_{\mathrm{ZV}}⟩$	$-0$	$.5731\left(42\right)$	$-0$	$.5770\left(29\right)$
$⟨F_P⟩$	$1$	$.1553\left(78\right)$	$1$	$.1643\left(35\right)$
$⟨F_{\mathrm{HFT}}+F_{\mathrm{ZV}}+F_P⟩$	$0$	$.5822\left(45\right)$	$0$	$.5860\left(35\right)$

Table 7. Table 7: Standard and tail-regression estimators of ⟨ A ⟩ delimited-⟨⟩ 𝐴 \langle A\rangle for the mixed DMC Hellmann-Feynman force on a carbon atom along the molecular axis of the C 2 molecule at r CC = 2 subscript 𝑟 CC 2 r_{\rm CC}=2 a.u., obtained from local force samples of various sizes M 𝑀 M . The “optimal” expansion orders n R subscript 𝑛 R n_{\rm R} and thresholds − log ⁡ q R subscript 𝑞 R -\log q_{\rm R} used for the tail-regression estimator in each case are also shown.

		Standard			TRE (unconstrained)					TRE (1 constraint)
$M$		$\overline{A}$			$n_R$	$-\log q_R$	$𝒜$			$n_R$	$-\log q_R$	$𝒜$
${10}^3$		$9$	$\left(35\right)$		$3$	$0.90$	$-14$	$\left(22\right)$		$3$	$1.10$	$-3$	$\left(11\right)$
${10}^4$		$-29$	$\left(28\right)$		$5$	$0.80$	$-9$	$.5\left(87\right)$		$4$	$0.90$	$-3$	$.3\left(33\right)$
${10}^5$		$4$	$.2\left(83\right)$		$4$	$0.90$	$-0$	$.3\left(28\right)$		$4$	$0.90$	$-0$	$.6\left(10\right)$
${10}^6$		$6$	$.1\left(38\right)$		$3$	$1.40$	$0$	$.9\left(11\right)$		$3$	$1.40$	$0$	$.53\left(45\right)$
${10}^7$		$-9$	$\left(12\right)$		$3$	$1.30$	$1$	$.08\left(34\right)$		$4$	$1.40$	$0$	$.61\left(18\right)$
${10}^8$		$2$	$.5\left(29\right)$		$3$	$1.40$	$0$	$.41\left(11\right)$		$4$	$1.40$	$0$	$.306\left(58\right)$

Table 8. Table 8: Standard and tail-regression estimators of ⟨ A ⟩ delimited-⟨⟩ 𝐴 \langle A\rangle for the mixed DMC local Hellmann-Feynman force and total force on a carbon atom along the molecular axis of the C 2 molecule at r CC = 2 subscript 𝑟 CC 2 r_{\rm CC}=2 a.u., obtained from 10 8 superscript 10 8 10^{8} sample points.

	Standard		TRE		TRE(1)
$⟨F_{\mathrm{HFT}}⟩$	$-2$	$.5\left(29\right)$	$-0$	$.41\left(11\right)$	$-0$	$.306\left(58\right)$
$⟨F_{\mathrm{HFT}}+F_{\mathrm{ZV}}+F_P⟩$	$0$	$.5810\left(42\right)$	$0$	$.5817\left(38\right)$

Equations

$$⟨A⟩=\frac{\int {\Psi }^{\ast }(\mathbf{R})\hat{A}\Psi (\mathbf{R})\mathrm{d}\mathbf{R}}{\int |\Psi (\mathbf{R}){|}^2\mathrm{d}\mathbf{R}},$$

$$P_{\mathbf{R}}(\mathbf{R})=\frac{|\Psi (\mathbf{R}){|}^2}{\int |\Psi (\mathbf{R}){|}^2\mathrm{d}\mathbf{R}},$$

$$⟨A⟩={\int }_{-\infty }^{\infty }{P}_A(A)A\mathrm{d}A,$$

$$P_A(A)={\int }_{\partial \Omega (A)}\frac{P_{\mathbf{R}}(\mathbf{R})}{|{\nabla }_{\mathbf{R}}A(\mathbf{R})|}{\mathrm{d}}^{dN-1}\mathbf{R},$$

$$\hat{H}(\mathbf{R})=-\frac{1}{2}\sum _i{\nabla }_i^2+\sum _{i,j>i}\frac{1}{r_{ij}}+\sum _{i,I}\frac{-Z_I}{r_{iI}}+\sum _{I,J>I}\frac{Z_I{Z}_J}{r_{IJ}},$$

$$P_E(E)=c_0|E-E_0{|}^{-4}+c_1|E-E_0{|}^{-5}+\dots ,$$

$$F(\mathbf{R})=F_{\mathrm{HFT}}(\mathbf{R})+F_{\mathrm{P}}(\mathbf{R}),$$

$$\begin{array}{rl}{F}_{\mathrm{HFT}}(\mathbf{R})& =-{\Psi }^{-1}(\mathbf{R})\frac{\partial \hat{H}(\mathbf{R})}{\partial x}\Psi (\mathbf{R})\\ & =\sum _i\frac{-Z_I{x}_{iI}}{r_{iI}^3}+\sum _{J\ne I}\frac{Z_I{Z}_J{x}_{IJ}}{r_{IJ}^3},\end{array}$$

$$F_{\mathrm{P}}(\mathbf{R})=-2{\Psi }^{-1}(\mathbf{R})\left[E(\mathbf{R})-⟨\hat{H}⟩\right]\frac{\partial \Psi (\mathbf{R})}{\partial x}.$$

$$F_{\mathrm{ZV}}(\mathbf{R})=-{\Psi }^{-1}(\mathbf{R})\left[\hat{H}(\mathbf{R})-E(\mathbf{R})\right]\frac{\partial \Psi (\mathbf{R})}{\partial x},$$

$$P_{F_{\mathrm{HFT}}}(F)=c_0|F-F_0{|}^{-5/2}+c_1|F-F_0{|}^{-3}+\dots ,$$

$$P_F(F)=c_0|F-F_0{|}^{-5/2}+c_1|F-F_0{|}^{-3}+\dots ,$$

$$⟨A⟩={\int }_{-\infty }^{\infty }{P}_A(A)A\mathrm{d}A,$$

$$\mathrm{Var}[A]={\sigma }^2={\int }_{-\infty }^{\infty }{P}_A(A){\left(A-⟨A⟩\right)}^2\mathrm{d}A.$$

$$\bar{A}=\frac{1}{M}\sum _{m=1}^M{A}_m,$$

$$S^2=\frac{{\sum }_{m=1}^M(A_m-\bar{A}{)}^2}{M-1}.$$

$$\mathrm{Var}[S^2]=\frac{1}{M}\left(m_4-\frac{M-3}{M-1}{\sigma }^4\right),$$

$$P_A(A)=\sum _n{c}_n{\left|A-A_0\right|}^{-\mu -n\Delta },$$

$$\begin{array}{rl}{P}_A(A)=& \sum _n{c}_n|A-A_{\mathrm{c}}-\epsilon {|}^{-\mu -n\Delta }\\ \approx & \sum _n{c}_n|A-A_{\mathrm{c}}{|}^{-\mu -n\Delta }\\ & +\epsilon \sum _n{c}_n(\mu +n\Delta )|A-A_{\mathrm{c}}{|}^{-\mu -n\Delta -1}\\ =& \sum _n{c}_n^{\prime }|A-A_{\mathrm{c}}{|}^{-\mu -n\Delta }=P_A^{\prime }(A),\end{array}$$

$$c_{n}^{\prime}=\left\{\begin{array}[]{l@{\,,\quad}l}c_{n}&n<\Delta^{-1}\\ c_{n}+\varepsilon c_{n-\Delta^{-1}}(\mu+n\Delta-1)&n\geq\Delta^{-1}\\ \end{array}\right.\;.$$

$$P_A(A)=\sum _{n=0}^{n_{\mathrm{R}}}{c}_n{\left|A-A_{\mathrm{c}}\right|}^{-\mu -n\Delta },A>A_{\mathrm{R}}.$$

$$⟨A⟩={\int }_{-\infty }^{A_{\mathrm{R}}}{P}_A(A)A\mathrm{d}A+{\int }_{A_{\mathrm{R}}}^{\infty }{P}_A(A)A\mathrm{d}A.$$

$$\mathcal{A}=\frac{1}{M}\sum _{m>M_{\mathrm{R}}}{A}^{(m)}+\sum _{n=0}^{n_{\mathrm{R}}}{c}_n{\int }_{A_{\mathrm{R}}}^{\infty }{\left|A-A_{\mathrm{c}}\right|}^{-\mu -n\Delta }A\mathrm{d}A.$$

$$\begin{array}{c}{\int }_{A_{\mathrm{R}}}^{\infty }|A-A_{\mathrm{c}}{|}^{-\mu -n\Delta }A\mathrm{d}A=\\ =\left[\frac{|A_{\mathrm{R}}-A_{\mathrm{c}}{|}^{-\mu -n\Delta +2}}{\mu +n\Delta -2}+A_{\mathrm{c}}\frac{|A_{\mathrm{R}}-A_{\mathrm{c}}{|}^{-\mu -n\Delta +1}}{\mu +n\Delta -1}\right].\end{array}$$

$$\mathcal{W}=\frac{M_{\mathrm{C}}}{M}+\sum _{n=0}^{n_{\mathrm{R}}}{c}_n{\int }_{A_{\mathrm{R}}}^{\infty }{\left|A-A_{\mathrm{c}}\right|}^{-\mu -n\Delta }\mathrm{d}A,$$

$$\begin{array}{rl}\mathcal{V}& =\frac{1}{M-1}\sum _{m>M_{\mathrm{R}}}(A^{(m)}-\mathcal{A}{)}^2\\ & +\sum _{n=0}^{n_{\mathrm{R}}}{c}_n{\int }_{A_{\mathrm{R}}}^{\infty }{\left|A-A_{\mathrm{c}}\right|}^{-\mu -n\Delta }(A-\mathcal{A}{)}^2\mathrm{d}A.\end{array}$$

$${\bar{F}}_A(A)={\int }_{A^{(m)}}^{\infty }{P}_A(A)dA\approx q_m,$$

$$\sum _{n=0}^{n_{\mathrm{R}}}\frac{c_n}{\mu +n\Delta -1}{\left|A^{(m)}-A_{\mathrm{c}}\right|}^{-\mu -n\Delta +1}=q_m,$$

$$q_m{\left|A^{(m)}-A_{\mathrm{c}}\right|}^{\mu -1}=\sum _{n=0}^{n_{\mathrm{R}}}\frac{c_n}{\mu +n\Delta -1}{\left|A^{(m)}-A_{\mathrm{c}}\right|}^{-n\Delta }.$$

$$\begin{array}{rl}{y}_m& =q_m{\left|A^{(m)}-A_{\mathrm{c}}\right|}^{\mu -1},\\ x_m& =|A^{(m)}-A_{\mathrm{c}}{|}^{-\Delta },\end{array}$$

Full text

I Introduction

Monte Carlo integration methods Kalos_Whitlock_montecarlo_1986 allow the evaluation of arbitrarily complicated high-dimensional integrals using random, discrete samples of the integrand. Besides the need to correct for serial correlation flyvbjerg_reblock_1989 ; Jonsson_reblock_2018 , the statistical analysis of these random samples is usually straightforward, and the final result of a Monte Carlo calculation is typically computed as a standard mean with an accompanying standard error, defining a confidence interval on the quantity of interest Stroock_probability_1993 . However there are problems for which the integrand diverges in such manner as to render these confidence intervals invalid.

This is the case in continuum quantum Monte Carlo (QMC) methods Ceperley_DMC_1980 ; Foulkes_QMC_2001 , a prominent family of tools for studying correlated many-body systems. Given a trial wave function $\Psi({\bf R})$, the variational Monte Carlo (VMC) method evaluates the expectation value of an observable $\hat{A}$ by accumulating its local value $A({\bf R})$ at random real-space configurations of the particles in the system, $\{\bf R\}$, distributed according to $\left|\Psi({\bf R})\right|^{2}$. The diffusion Monte Carlo (DMC) method samples the lowest-energy wave function $\Phi({\bf R})$ with the same nodal structure as $\Psi({\bf R})$ by stochastic projection according to the imaginary-time Schrödinger equation, and yields more accurate estimates of observables than VMC. The stochastic integration employed by these methods allows using trial wave functions that are not analytically integrable, providing extraordinary flexibility and compactness in the description of many-body correlations LopezRios_Jastrow_2012 ; LopezRios_backflow_2006 ; Bajdich_Pfaffians_2006 . The VMC and DMC methods are routinely used to solve electronic structure and quantum chemistry problems Ceperley_DMC_1980 ; Seth_atoms_2011 ; Drummond_hydrogen_2015 .

The mismatch between the nodes of the trial wave function, $\{{\bf R}:~{}\Psi({\bf R})=0\}$, and those of the true ground-state wave function of the system are the main source of outliers in the local energy distribution sampled in QMC, resulting in heavy tails Panharipande_helium_1986 ; Trail_htail_2008 that preclude the evaluation of meaningful confidence intervals on the estimated variance of the local energy. The local atomic force, comprising a Hellmann-Feynman force Feynman_forces_1939 and a Pulay term Pulay_forces_1969 ; Casalegno_forces_2003 , has heavy tails arising both from the divergence of the electron-nucleus potential energy in all-electron systems Badinski_forces_2010 and from the nodal error, which prevent the evaluation of meaningful confidence intervals on the estimated expectation value of the force.

Various methods have been proposed to circumvent the statistical hurdles in the evaluation of atomic forces in QMC. Modified estimators of the force satisfying a zero-variance principle have been proposed Assaraf_zero-variance_1999 ; Assaraf_forces_2001 ; Assaraf_forces_2003 ; Per_forces_2008 ; Badinski_zero-variance_2008 that substantially reduce the magnitude of the heavy tails in the local force distribution. Force estimation methods based on the use of pseudopotentials Badinski_pp_forces_2007 ; Badinski_pp_forces_2008 can eliminate the problematic behavior of the Hellmann-Feynman force Badinski_forces_2010 , as does the fitting approach proposed by Chiesa et al. Chiesa_force_density_2005 . However, some of these methods involve approximations, and none of them addresses the heavy tails in the local Pulay force distribution, and therefore the total force remains affected by an infinite-variance problem. The reweighted approach proposed by Attaccalite and Sorella Attaccalite_reweight_2008 does prevent the Pulay force from diverging, but it involves modifying the sampling distribution and is therefore not applicable in DMC.

Tail-index estimation methods Hill_TIE_1975 ; Beirlant_TIE_1996 ; Kratz_qq_1996 ; Baek_TIE_2010 allow the exponent governing the asymptotic behavior of heavy-tailed probability distributions to be estimated from statistical samples. This prior work, combined with knowledge of the exact asymptotic form of the tails of the local energy Trail_htail_2008 and local force Badinski_forces_2010 distributions, provides us with the foundation to develop the tail-regression estimator (TRE) for heavy-tailed distributions. We demonstrate the application of this technique to VMC and DMC data, for which we are able to obtain forces and local energy variances not affected by the infinite-variance problem.

The rest of this paper is structured as follows. We review the properties of the heavy tails of the local energy and local force distributions in section II. In section III we discuss the standard method for estimating an expectation value from a statistical sample, and we propose our tail-regression estimator in section IV. The application of our methodology is illustrated in section V using model distributions of known statistical properties. Finally, we present the results of applying our method to the VMC energy of a homogeneous electron gas and the VMC and DMC atomic force on an atom in the all-electron C2 molecule in sections VI and VII, and our conclusions are stated in section VIII.

II Heavy tails in quantum Monte Carlo

We explore the formal definition of expectation values in QMC to allow the characterization of the resulting heavy-tailed distributions. For simplicity, we restrict our analysis to the VMC method until Section VII, in which we discuss the application of our methodology to DMC data. Note that we ignore serial correlation in this work.

Given a trial wave function $\Psi({\bf R})$, a VMC calculation evaluates the expectation value

[TABLE]

by generating electronic configurations $\{\bf R\}$ according to the distribution

[TABLE]

and evaluating the local values $A({\bf R})=\Psi^{-1}({\bf R})\hat{A}\Psi({\bf R})$ of observable $\hat{A}$. The VMC expectation value can be recast as

[TABLE]

where

[TABLE]

where $d$ is the dimensionality of the system, $N$ is the number of electrons, and $\partial\Omega(A)$ is the $(dN-1)$-dimensional region of configuration space where $A({\bf R})=A$.

The local value of some important observables diverge at certain configurations, and it is often possible to characterize the asymptotic behavior of $P_{A}(A)$ from knowledge of the analytical form of $\Psi({\bf R})$ near to these configurations. We summarize the relevant properties of the local energy and the local atomic force below.

II.1 Local energy

Consider the Hamiltonian of a molecular system in Hartree atomic units ($\hbar=m_{\rm e}=|e|=4\pi\epsilon_{0}=1$),

[TABLE]

where $r_{ij}$ is the distance between the $i$th and $j$th electrons, $r_{iI}$ is the distance between the $i$th electron and the $I$th nucleus, $r_{IJ}$ is the distance between the $I$th and $J$th nuclei, and $Z_{I}$ is the atomic number of the $I$th nucleus.

The situations in which the local energy $E({\bf R})=\Psi^{-1}({\bf R})\hat{H}({\bf R})\Psi({\bf R})$ diverges were classified in detail by Trail Trail_htail_2008 . The divergence of the Coulomb potential at electron-nucleus and electron-electron coalescence points, $r_{iI}\to 0$ and $r_{ij}\to 0$ respectively, can be neutralized by constraining the trial wave functions to obey the Kato cusp conditions Kato_cusp_1957 under which the kinetic energy exactly cancels the potential energy at two-body coalescence points. The remaining divergences of the local energy arise when $\Psi({\bf R})\to 0$ but $\Psi$ is not locally identical to an eigenstate of $\hat{H}$, since $\hat{H}({\bf R})\Psi({\bf R})$ can be finite where $\Psi({\bf R})$ is zero.

Mismatches between the nodes of $\Psi({\bf R})$ and $\hat{H}({\bf R})\Psi({\bf R})$ are responsible for the asymptotic behavior Trail_htail_2008

[TABLE]

when $|E|\to\infty$, where $E_{0}$ is the exact ground state energy and $\{c_{i}\}$ are unknown coefficients. The coefficients of even powers of $|E-E_{0}|$ in the left ($E\to-\infty$) and right ($E\to+\infty$) tails of $P_{E}(E)$ have equal coefficients, $c_{2n}^{\rm L}=c_{2n}^{\rm R}$, while those of odd powers are of the same magnitude but of opposite signs, $c_{2n+1}^{\rm L}=-c_{2n+1}^{\rm R}$.

The expectation value of the energy itself can be evaluated with standard estimators without problems, but these yield unreliable confidence intervals on the variance of the local energy. The variance of the local energy is an important quantity since it is directly related to the quality of the trial wave function, and is routinely used in wave function optimization Umrigar_varmin_1988 , as well as in the “variance extrapolation” technique that attempts to estimate the zero-variance (exact wave function) limit of expectation values Holzmann_ueg_2003 . We discuss the specific issues with the variance of the local energy in section III.

II.2 Local force

The force exerted by the electrons and the other nuclei on the $I$th nucleus of a system along the $x$ direction is $\langle\hat{F}_{x,I}\rangle=-{\rm d}\langle\hat{H}\rangle/{\rm d}x_{I}$, where $x_{I}$ is the $x$ Cartesian coordinate of the $I$th nucleus. Dropping the $I$ and $x$ labels, the local force can be expressed as

[TABLE]

where the Hellmann-Feynman force is

[TABLE]

and the Pulay force is

[TABLE]

Optionally, a variance-reduction term which satisfies a zero-variance principle Assaraf_forces_2003 ,

[TABLE]

can be added to Eq. 7. This term does not alter the expectation value of the force but reduces the extent of the fluctuations of the local Hellmann-Feynman force.

The local Hellmann-Feynman force diverges at electron-nucleus coalescence points, and its distribution exhibits a power-law tail of the form

[TABLE]

when $|F|\to\infty$, where $F_{0}$ is a constant and $\{c_{i}\}$ are unknown coefficients. The coefficients of the leading-order term on the left and right tails are equal, $c_{0}^{\rm L}=c_{0}^{\rm R}$, and the rest are asymmetric.

If the wave function satisfies the electron-nucleus Kato cusp conditions Ma_cusp_2005 the zero-variance term exactly cancels this divergence Per_forces_2008 , but, like for the local energy, the mismatch between the nodes of $\Psi(\bf R)$ and $\hat{H}({\bf R})\Psi({\bf R})$ is responsible for the heavy tails in the distribution of the local values of the Pulay and zero-variance terms, and hence of the total force, which satisfies Badinski_forces_2010

[TABLE]

when $|F|\to\infty$, where $F_{0}$ is a constant and $\{c_{i}\}$ are unknown coefficients exhibiting no symmetry.

A somewhat different scenario arises if nondivergent pseudopotentials are used in place of the electron-nucleus Coulomb potential Trail_pp1_2015 ; Trail_pp2_2017 and the local force estimation is consequently adjusted Badinski_pp_forces_2007 ; Badinski_pp_forces_2008 . In this case the local Hellmann-Feynman force exhibits less problematic heavy tails of leading order $|F-F_{0}|^{-4}$ Badinski_forces_2010 , but since the Pulay term is unaffected by the use of pseudopotentials the local total force remains of the form of Eq. 12.

III Standard estimation of an expectation value

The expectation value of an observable whose local value $A$ is distributed according to $P_{A}(A)$ is

[TABLE]

and the variance of $A$ is the expectation value of $\left(A-\langle A\rangle\right)^{2}$,

[TABLE]

The integrals in Eqs. 13 and 14 must be nondivergent for the expectation value and variance of $A$ to be well-defined. Therefore, probability distributions with asymptotic behavior $P_{A}(A)\sim|A|^{-\mu}$ as $|A|\to\infty$ have no well-defined expectation value or variance for $\mu\leq 2$, and have a well-defined expectation value but no well-defined variance for $2<\mu\leq 3$. Note that a function with $\mu\leq 1$ is not a valid probability distribution as it cannot be normalized.

Let $\{A_{m}\}_{m=1}^{M}$ be a sample of $M$ independent random variables identically distributed according to $P_{A}(A)$. The standard estimator for Eq. 13 is the sample mean,

[TABLE]

and the standard estimator for Eq. 14 is the sample variance,

[TABLE]

The uncertainty on $\bar{A}$ is the standard error $\sigma_{\bar{A}}=S/\sqrt{M}$. This poses a problem for distributions of leading-order exponent $2<\mu\leq 3$: despite having a well-defined expectation value according to Eq. 13, its standard estimator has a divergent uncertainty because it is defined in terms of the divergent variance of Eq. 14.

In this regime $P_{A}(A)$ does not satisfy the conditions of the central limit theorem that would guarantee the asymptotic normality of confidence intervals built from the standard mean and standard error Stroock_probability_1993 . Instead, $P_{A}(A)$ satisfies the law of large numbers, which states that the standard mean does converge to the expectation value at infinite sample size but confidence intervals cannot be constructed using the standard error as finite sample sizes.

A similar issue affects the estimator of the variance itself. Even though the variance is well-defined for $\mu>3$, the variance on the estimator of the variance is

[TABLE]

where $m_{4}$ is the fourth-order central moment of $P_{A}(A)$, which diverges for distributions of leading-order exponent $\mu\leq 5$, leading to a divergent uncertainty on the $S^{2}$ estimator of the variance.

IV Tail-regression estimator

As outlined in the previous section, the uncertainty on the standard estimator of a moment of a probability distribution involves the estimator of a higher-order moment, which might be divergent even though the moment of interest is well-defined. We therefore propose an alternative estimator of the moment of a heavy-tailed probability distribution that exploits knowledge of its analytical asymptotic form to yield well-defined confidence intervals whenever the moment itself is well defined.

Without loss of generality we focus on distributions with a right heavy tail; the extension of our analysis to distributions with left and both left and right heavy tails is straightforward. In particular, we consider a probability distribution exhibiting a right tail of asymptotic form

[TABLE]

when $A\to\infty$, where the leading-order exponent $\mu$ and exponent increment $\Delta>0$ are assumed to be known analytically, as is the case for local energies, with $\mu=4$ and $\Delta=1$ Trail_htail_2008 , and for local forces, with $\mu=5/2$ and $\Delta=1/2$ Badinski_thesis_2008 . The specific value of $\Delta$ is not critical for the correct description of the asymptote by Eq. 18, and can be assumed to be unity, but the accuracy of a truncated expansion strongly depends on $\Delta$. In other words, the bias incurred by choosing a suboptimal value of $\Delta$ can be made arbitrarily small by using a larger expansion. In Eq. 18 $\{c_{n}\}$ are unknown coefficients, and $A_{0}$ is an unknown parameter which is assumed to lie close to the “center” of the distribution.

IV.1 Validity of asymptote with approximate $A_{0}$

First, we address the fact that $A_{0}$ is an unknown nonlinear parameter in Eq. 18 and will have to be approximated. Let $A_{\rm c}$ be an approximation to $A_{0}$ such that $A_{0}=A_{\rm c}+\varepsilon$, where $\varepsilon$ is a small error. Assuming for simplicity that $\Delta^{-1}$ is an integer and expanding to first order in $\varepsilon$ we find

[TABLE]

where

[TABLE]

The asymptotic expression $P_{A}^{\prime}(A)$ has the same form as $P_{A}(A)$, albeit with modified coefficients $c_{n}^{\prime}$ for $n\geq\Delta^{-1}$. We shall therefore proceed with the derivation of our estimator using $P_{A}^{\prime}(A)$ as the asymptote, dropping the primes from the notation for clarity. This effectively amounts to replacing $A_{0}$ with $A_{\rm c}$, which in practice we set to the sample median.

IV.2 Estimator

In order to develop our estimator, we start by assuming that there exists a threshold $A_{\rm R}$ such that for $A>A_{\rm R}$ the probability distribution is accurately represented by an expansion of order $n_{\rm R}$,

[TABLE]

The integral of Eq. 13 can be partitioned at $A_{\rm R}$ into central and right-tail contributions,

[TABLE]

Let $\{A_{m}\}$ be a sample of $M$ independent random variables identically distributed according to $P_{A}(A)$, $\{A^{(m)}\}$ the corresponding order statistics, i.e., the re-indexed version of $\{A_{m}\}$ such that $A^{(1)}>A^{(2)}>\ldots>A^{(M)}$, $M_{\rm C}$ the number of data in the central region $A<A_{\rm R}$, and $M_{\rm R}=M-M_{\rm C}$ the number of data in the tail. We define the tail-regression estimator of $\langle A\rangle$ as

[TABLE]

The integrals in Eq. 23 are nondivergent for $\mu>2$ and can be evaluated analytically,

[TABLE]

The parameters $\{c_{n}\}$ in Eq. 23 can be obtained by regression of $\{A^{(m)}\}_{m=1}^{M_{\rm R}}$ to Eq. 21. Since $\cal A$ is linear in $\{c_{n}\}$, the distribution of $\cal A$ follows that of $\{c_{n}\}$. This implies that if the regression coefficients are asymptotically normally distributed, ${\cal A}$ will also be asymptotically normally distributed. We will address the distribution of regression coefficients in section IV.4. We regard $A_{\rm R}$, and $n_{\rm R}$ as external parameters that we deal with separately, see section IV.5, and do not contribute to the uncertainty on ${\cal A}$.

Analogously, we define the tail-regression estimators of the norm,

[TABLE]

and of the variance of the distribution,

[TABLE]

These integrals can likewise be evaluated analytically. Finally, we note that the threshold $A_{\rm R}$ must lie on the midpoint between two adjacent sample points, $\frac{1}{2}\left(A^{(M_{\rm R})}+A^{(M_{\rm R}+1)}\right)$, to ensure that the central contributions have the correct weight in Eqs. 23, 25, and 26.

We use the bootstrap method Efron_bootstrap_1986 to compute the uncertainty on ${\cal A}$. We generate $n_{\rm bs}$ resamples of $\{A_{m}\}$ with replacement, that is, the $i$th resample $\{A_{m}^{[i]}\}$ contains $M$ elements from $\{A_{m}\}$ drawn at random and uniformly, allowing repetitions. For each resample we evaluate ${\cal A}^{[i]}$, and we evaluate the uncertainty on ${\cal A}$ as the standard deviation of the values of $\{{\cal A}^{[i]}\}$. Estimates of other statistical parameters arising from analysis of $\{A_{m}\}$, including $\cal W$ and $\cal V$, are also obtained in this process. In our applications of the tail-regression estimator we use $n_{\rm bs}=4096$ bootstrap resamples, which provide a $1.1\%$ uncertainty on the estimated uncertainties assuming normality. The computational cost of this approach is proportional to $M\times n_{\rm bs}$.

IV.3 Tail form in $yx$ scale

The framework for our tail-regression procedure is inspired by regression-based tail-index estimation methods Hill_TIE_1975 ; Beirlant_TIE_1996 , discussed in Appendix A. First, note that the complementary cumulative distribution function $\bar{F}_{A}(A)$ associated with $P_{A}(A)$ should be approximately equal to the sample quantiles $q_{m}$,

[TABLE]

for which we use the symmetric form $q_{m}=\frac{m-1/2}{M}$. Substituting Eq. 21 into Eq. 27 yields

[TABLE]

which we rearrange as

[TABLE]

We define

[TABLE]

which we refer to as “$yx$ scale”, under which Eq. 29 reads

[TABLE]

that is, $y$ is simply a polynomial in $x$. By construction, $y$ must be positive and tend to a finite value as $x\to 0$, $y(0)=c_{0}/(\mu-1)$, and the $n$th derivative of $y(x)$ at $x=0$ is likewise proportional to $c_{n}$.

It is useful to inspect the basic properties of the $yx$ scale we have introduced. For illustration purposes, let

[TABLE]

which for $\mu>2$ is a normalized probability distribution whose expectation value is zero and has the asymptote $|A|^{-\mu}$ as $|A|\to\infty$. In Fig. 1 we show a $yx$-scale plot of the left tail of $M$ independent random numbers identically distributed according to $H_{4}(A)$ at different sample sizes $M$, assuming $\Delta=1$. The exact value $y(0)=0.1501$ is shown as a short-dashed line in each panel.

The first quantile $q_{1}=1/(2M)$, corresponding to the largest value of $A$ in the sample, gets smaller as $M$ increases, and, as follows from Eq. 30, the extreme point $(x_{1},y_{1})$ satisfies $y_{1}=q_{1}x_{1}^{(1-\mu)/\Delta}$. This curve determines how far left the plot extends, as shown by the long-dashed line in each of the panels of Fig. 1. With increasing $M$ the plot is populated from right to left, gradually producing a better resolved curve near $x=0$.

In Fig. 2 we show $yx$ plots of the $M=10^{6}$ sample used in Fig. 1 in which the $yx$ scale is defined using the exponents $\mu^{\prime}=3$, $4$, and $5$.

If the true exponent is underestimated, $\mu^{\prime}<\mu$, the curve goes to zero as $x\to 0$, while overestimating the exponent, $\mu^{\prime}>\mu$, yields a diverging curve at $x\to 0$. Attempting to use $yx$-scale plots to estimate $\mu$ is rudimentary at best, since each possible value must be tested explicitly, and care must be taken in interpreting the large statistical fluctuations of $y$ at small $x$. Tail-index estimation methods Beirlant_TIE_1996 ; Kratz_qq_1996 offer a robust alternative; the regression method of Ref. Beirlant_TIE_1996 produces a reasonable value of $\mu=3.94(2)$ from the data plotted in Fig. 2. Note that the tail-regression estimator relies on analytical knowledge of $\mu$ rather than on its estimation from the sample.

IV.4 Tail regression and weights

We now focus on the regression procedure to ensure a reliable estimation of the fit parameters. It is apparent from Fig. 1 that the distribution of $y$ values at small $x$ about the exact $y(0)$ is skewed towards large $y$, so adequate least-squares weights are needed to ensure the faithfulness of the resulting fit.

Let $r_{m}=w_{m}[y_{m}-y(x_{m})]$ be the least-square residuals, where $w_{m}$ is the least-squares weight applied to the $m$th point, and $\chi^{2}=\frac{1}{M-n_{\rm R}-1}\sum_{m=1}^{M}r_{m}^{2}$ be the least-squares function. Minimizing $\chi^{2}$ with respect to the fit parameters $\{c_{n}\}$ yields the linear system of equations $T{\bf c}={\bf b}$, where

[TABLE]

and indices $p$ and $q$ run between [math] and $n_{\rm R}$. The parameter vector is thus

[TABLE]

where ${\rm adj}(T)$ is the adjugate matrix of $T$ (i.e., the transpose of its cofactor matrix) and $\det(T)$ is its determinant.

The parameters $({\bf c})_{p}$ will be asymptotically normally distributed if $(T)_{pq}$, $({\bf b})_{p}$, and $\det(T)$ are themselves asymptotically normally distributed and $\det(T)$ is nonzero, since Eq. 34 involves sums and multiplications, which preserve asymptotic normality, and a division, which also preserves asymptotic normality if the denominator is strictly nonzero.

We consider weights of the form

[TABLE]

where $\gamma$ is a positive constant. Since $x_{m}$ is a negative power of $|A^{(m)}-A_{\rm c}|$, all elements of $T$ and $\det(T)$ are asymptotically normally distributed, but the elements of $\bf b$ need not be. We focus on the $m=1$ contribution to $({\bf b})_{0}$, which is the least likely to exhibit asymptotic normality,

[TABLE]

Noting that $w_{1}=q_{1}^{\gamma}y_{1}^{-\gamma}$ and that the exact asymptotic value of $y_{1}$ as $M\to\infty$ is by construction $y(0)=c_{0}/(\mu-1)$, the asymptotic limit of $w_{1}y_{1}$ is

[TABLE]

where $q_{1}=1/(2M)$ encodes the sample-size dependence.

We now investigate the distribution of values of $w_{1}y_{1}$ about $\xi$. The probability that $|A^{(1)}-A_{\rm c}|$ is bounded from above by $\alpha$ is the probability that the $M$ points in the sample are bounded from above by $A_{\rm c}+\alpha$, that is,

[TABLE]

to leading order for large $\alpha$. Differentiating with respect to $\alpha$ yields the probability distribution of the extreme value,

[TABLE]

and by a change of variable,

[TABLE]

assuming $\gamma\neq 1$. At large $\beta$,

[TABLE]

which is a power-law tail for $\gamma<1$, yielding an undefined expectation value for $\gamma\leq 1/2$. Unweighted fits ($\gamma=0$) are therefore numerically ill-conditioned since the residuals are themselves heavy tailed $\sim\beta^{-2}$ regardless of the value of $\mu$. Fit weights with $\gamma>1/2$ must therefore be used.

In Fig. 3 we plot $P_{w_{1}y_{1}/\xi}(\beta)$ for various values of $\gamma$. Note that the curve for $\gamma=0$ describes the distribution of values of $y$ in Fig. 1 relative to the analytical value of $y(0)$ along the constant-$q$ lines. The plots for $\gamma=0.8$ and $1.2$ reveal that, while not ill-conditioned, the significant skewness of the distributions makes the first moment of $P_{w_{1}y_{1}/\xi}(\beta)$ differ significantly from the asymptotic expectation value of $1$ in these cases. For $\gamma=1$, $w_{1}y_{1}=q_{1}^{\gamma}$ is $y_{1}$-independent, and the distribution is therefore a delta function peaked at the asymptotic expectation value, as shown in Fig. 3. We therefore use $\gamma=1$ for our fit weights.

We empirically find it advantageous to include a $q_{m}$-dependent factor in the fit weights so as to ensure the continuity of the fit to the probability distribution near $A_{\rm R}$. We choose this factor to be the weights corresponding to the formulation of the Hill estimator of the first-order tail index Hill_TIE_1975 as a regression estimator Beirlant_TIE_1996 , see Eq. 50 in Appendix A. Therefore, our full fit weights are

[TABLE]

IV.5 Selecting $n_{\rm R}$ and $A_{\rm R}$

The tail-regression estimator depends parametrically on the expansion order $n_{\rm R}$ and threshold $A_{\rm R}$. In our tests we try several thresholds and converge the fit with respect to the expansion order at each of them. We then choose the value of $A_{\rm R}$ which minimizes the uncertainty on either $\cal V$, if well-defined, or on $\cal A$.

We choose the expansion order heuristically by finding plateaus in $\cal W$, $\cal A$, $\cal V$, and $\chi^{2}$ as a function of $n_{\rm R}$, and selecting the smallest expansion order at which these four functions have converged. To ensure correctness, we further require that ${\cal W}\approx 1$ and $y(x)>0$ within the fit range, and we restrict $n_{\rm R}\geq 1/\Delta$ in order to “absorb” the error incurred by approximating $A_{0}$ by $A_{\rm c}$, as explained in Section IV.1.

Note that we do not set $A_{\rm R}$ directly, but instead set $q_{\rm R}=(M_{\rm R}-1/2)/M$, as keeping the number of sample points in each partition fixed across bootstrap resamples eliminates the variation of the central contribution to $\cal W$, which is statistically advantageous. We choose our values of $q_{\rm R}$ using a grid of equally-spaced values of $-\log q_{\rm R}$. We pick $n_{\rm R}$ and $A_{\rm R}$ using $n_{\rm bs}=256$, and evaluate the final result separately with $n_{\rm bs}=4096$ to avoid selection bias. We illustrate the procedure for choosing $n_{\rm R}$ and $A_{\rm R}$ in section V.1.

IV.6 Two-tailed distributions, symmetry, and constraints

The tail-regression estimator described so far can be modified trivially for distributions with left and right heavy tails. For simplicity we use the same expansion orders and thresholds on both tails, $n_{\rm R}=n_{\rm L}$ and $M_{\rm R}=M_{\rm L}$.

In some important cases, including the local energy and local atomic force in VMC Trail_htail_2008 ; Badinski_forces_2010 , the leading-order coefficients of the left and right tails, $c_{0}^{\rm L}$ and $c_{0}^{\rm R}$, are equal. This can be exploited by unifying the regression step for both tails and imposing the constraint $c_{0}^{\rm L}=c_{0}^{\rm R}=c_{0}$. The leading order contribution to $\cal A$ from the tails is

[TABLE]

The exact cancellation of part of the left- and right-tail contributions to $\cal A$ should provide a substantial reduction to its uncertainty. The effect on the uncertainty on $\cal V$ of enforcing symmetry can be expected to be marginal since both tails contribute positively in this case.

As implied by Eq. 20, due to the approximation $A_{0}\approx A_{\rm c}$ constraints must not be applied to parameters $c_{n}$ with $n\geq 1/\Delta$. For example, even if a distribution with $\Delta=1$ is known to analytically satisfy $c_{1}^{\rm L}=c_{1}^{\rm R}$, the values of the $c_{1}$ parameters on each tail must be allowed to differ to account for the error in $A_{\rm c}$. In the tests carried out in our present work we use at most one constraint, and we use the labels “TRE” and “TRE(1)” in the plots in sections V, VI, and VII to distinguish the unconstrained and constrained estimators.

The use of constraints allows for the interesting possibility of estimating the expectation value of distributions with $1<\mu\leq 2$ for which $\langle A\rangle$ is formally divergent. In this case we redefine the expectation value as the Cauchy principal value of the integral in Eq. 13 with respect to $A_{0}$,

[TABLE]

ensuring that the divergent leading-order contributions cancel out due to symmetry. We present an example of this in section V.3.

V Application to model distributions

In this section we apply the tail-regression estimator to synthetic data. We construct the seed model distributions as linear combinations of $H_{\mu}(A)$, defined in Eq. 32, to study various cases of interest.

V.1 Distribution with undefined fourth moment

Distributions with $3<\mu\leq 5$ have convergent standard estimators for the expectation value and variance, but the uncertainty on the standard estimator of the variance is divergent. To exemplify this case we choose to analyze $P_{A}(A)=\frac{1}{2}H_{3.1}(A)+\frac{1}{2}H_{4.1}(A)$, which has a leading-order tail exponent of $\mu=3.1$, close to the lower limit of $3$, and $\Delta=1$. The analytical variance of this distribution is $\sigma^{2}=4.6586$, and $P_{A}(A)$ satisfies the analytical limit $y\sim 0.0997+0.0730x$ as $x\to 0$.

The regression of the tails of this model distribution is demonstrated in Fig. 4 using a sample of $10^{6}$ random variables. The top panel shows the probability distribution estimated by convolving the data with a variable-width Gaussian kernel, and the lower panels show $yx$ plots of the data. The tail fits are shown in the three panels as thick lines, for which we use $n_{\rm R}=3$ and $-\log q_{\rm R}=2.25$, and the constraint $c_{0}^{\rm L}=c_{0}^{\rm R}$ has been imposed.

The process of selecting $n_{\rm R}$ at fixed threshold is illustrated in Fig. 5, where $\cal W$, $\cal A$, $\cal V$, and $\chi^{2}$ are plotted as a function of $n_{\rm R}$. These functions converge relatively quickly with $n_{\rm R}$, but for large expansion orders overfitting becomes an issue. This can be easily spotted by the significant jump in the uncertainty on the estimators, which in Fig. 5 occurs at $n_{\rm R}=7$. In this case we consider $n_{\rm R}=3$ to be the “optimal” expansion order.

In Fig. 6 we plot the uncertainty on $\cal V$ as a function of the threshold for each considered expansion order. The uncertainty on $\cal V$ is largely monotonic in $n_{\rm R}$ at fixed threshold and in $-\log q_{\rm R}$ at fixed expansion order. Figure 6 shows that the “optimal” values of the expansion order and threshold are $n_{\rm R}=3$ and $-\log q_{\rm R}=2.25$, as used in Fig. 4.

In Table 1 we report the results obtained from the standard and tail-regression estimators of $\langle A\rangle$ and $\sigma^{2}$ for sample sizes ranging from $10^{3}$ to $10^{8}$, which we plot in Figs. 7 and 8.

For the tail-regression estimator we report results both without the use of constraints and imposing $c_{0}^{\rm L}=c_{0}^{\rm R}$. The “optimal” values of $n_{\rm R}$ and $-\log q_{\rm R}$ are not significantly sample-size dependent; we find that $n_{\rm R}$ increases when $-\log q_{\rm R}$ decreases and vice versa, as could be expected.

The estimators of $\langle A\rangle$ are within uncertainty of the exact value of zero at all sample sizes, but $\cal A$ has an uncertainty about $25\%$ smaller than $\bar{A}$. The uncertainty on the standard estimator of the variance $S^{2}$ is nonmonotonic, as expected, and in this case confidence interval sizes are severely underestimated, causing the false impression that $S^{2}$ converges to an incorrect value with increasing $M$. By contrast, $\cal V$ remains within uncertainty of $\sigma^{2}$ at all $M$, and its uncertainty decreases monotonically with sample size. The unconstrained and constrained estimators give indistinguishable results in this case, and the uncertainty on both seems to asymptotically decay as $M^{-1/2}$.

V.2 Distribution with undefined second moment

Distributions with $2<\mu\leq 3$ have a divergent variance, and the standard estimator of $\langle A\rangle$ has an undefined uncertainty. We exemplify this case with model distribution $P_{A}(A)=\frac{1}{2}H_{2.1}(A)+\frac{1}{2}H_{3.1}(A)$, which has $\mu=2.1$, close to the lower limit of $2$, and $\Delta=1$.

The standard and tail-regression estimators of $\langle A\rangle$ are given in Table 2 and plotted in Fig. 9 as a function of sample size $M$.

As expected, the standard estimator $\bar{A}$ hovers around the exact value $\langle A\rangle=0$ but its uncertainty does not decrease uniformly with sample size. The tail-regression estimator provides a monotonically decreasing uncertainty with an approximate asymptotic decay proportional to $M^{-1/2}$, and imposing the constraint $c_{0}^{\rm L}=c_{0}^{\rm R}$ yields an order of magnitude smaller uncertainties than the unconstrained estimator does.

V.3 Symmetric distribution with undefined first moment

Distributions with $1<\mu\leq 2$ have a divergent expectation value, and the standard estimator of $\langle A\rangle$ is undefined. However if the tails of the distribution are symmetric to leading order it is possible to redefine $\langle A\rangle$ as a Cauchy principal value which can be estimated, see Eq. 44. We exemplify this case with model distribution $P_{A}(A)=\frac{1}{2}H_{1.1}(A)+\frac{1}{2}H_{2.1}(A)$, which has $\mu=1.1$, close to the lower limit of $1$, and $\Delta=1$.

Results using the constrained tail-regression estimator are given in Table 3 and plotted in Fig. 10.

We find that $\cal A$ is within uncertainty of the exact value of zero at all sample sizes, and again the uncertainty on $\cal A$ appears to be asymptotically proportional to $M^{-1/2}$. In Table 3 we also give the order of magnitude of the computed sample mean to highlight that this example is absolutely intractable with the standard estimator.

VI Application to variational quantum Monte Carlo data

In this section we explore the performance of the tail-regression estimator on data obtained from VMC calculations. Local values of observables generated using the VMC method are usually affected by serial correlation due to the use of the Metropolis algorithm to sample configuration space. In our VMC calculations we perform up to $400$ Metropolis steps between consecutive evaluations of the local values of the target observables so that these can be considered to be independent random variables, and we have verified that the resulting datasets exhibit negligible serial correlation.

VI.1 The energy of the homogeneous electron gas

The homogeneous electron gas is an ideal test bed for methodological developments in QMC. We perform VMC calculations on the paramagnetic 54-electron gas in a cubic simulation cell at density $r_{\rm s}=1$ using the Slater-Jastrow (SJ) wave function, consisting of the product of up- and down-spin Slater determinants of the plane-waves with the smallest momenta compatible with the periodicity of the simulation cell multiplied by a Jastrow correlation factor, $\Psi_{\rm SJ}({\bf R})={\rm e}^{J({\bf R})}D_{\uparrow}({\bf R}_{\uparrow})D_{\downarrow}({\bf R}_{\downarrow})$. Our Jastrow factor consists of an isotropic electron-electron term of the Drummond-Towler-Needs form Drummond_Jastrow_2004 ; LopezRios_Jastrow_2012 . We use the casino code casino_reference to generate samples of local energies whose distribution, as detailed in section II.1, has left and right heavy tails of principal exponent $\mu=4$, $\Delta=1$, and equal left- and right-tail leading-order coefficients Trail_htail_2008 . As explained in section IV.6, constraints involving $c_{n}$ with $n\geq 1$ cannot be applied since $A_{0}$ is being approximated by $A_{\rm c}$.

In Fig. 12 we plot the probability distribution estimated from $10^{6}$ local energies, $yx$ plots of the tails, and the corresponding tail fits.

The standard and tail-regression estimators of $\langle A\rangle$ and $\sigma^{2}$, given in Table 4 and plotted in Fig. 11 as a function of sample size, are in good agreement with each other.

The $\cal A$ estimator offers no advantage over $\bar{A}$ in this case, but the uncertainty on $S^{2}$ exhibits nonmonotonicity as a function of $M$, while that in $\cal V$ is monotonic, significantly smoother, and up to 45% smaller. Note that even though the nominal standard confidence interval on $S^{2}$ only shows minor signs of ill behavior in this example, it is formally undefined, while the tail-regression estimator produces valid confidence intervals. This is of potential practical importance in wave function optimization and variance extrapolation.

VI.2 The atomic force in the C2 molecule

We turn our attention to the atomic force in the all-electron carbon dimer. The C2 molecule is of particular interest due to its strong multi-reference character that makes the single-determinantal wave function incur a large nodal error, which ought to provide relatively strong heavy tails in the local Pulay and zero-variance force distributions.

We generate Hartree-Fock orbitals for the all-electron carbon dimer at an off-equilibrium (compressed) bond length of $r_{\rm CC}=2.0$ a.u. (the experimental equilibrium bond length of C2 is $2.3481$ a.u. Toulouse_Umrigar_2008 ) using the relatively modest cc-pvdz basis set Dunning_basis_1992 ; Kendall_basis_1992 with the molpro code molpro . We combine these orbitals, modified to satisfy the Kato cusp conditions at electron-nucleus coalescence points Ma_cusp_2005 , with a Jastrow factor containing electron-electron, electron-nucleus, and electron-electron-nucleus terms of the Drummond-Towler-Needs form Drummond_Jastrow_2004 ; LopezRios_Jastrow_2012 , to form the trial wave function for VMC. Throughout the VMC run, performed with the casino code casino_reference , we collect local values of the components of the force on one of the carbon atoms along the molecular axis in the direction away from the other atom.

We first focus on the Pulay force which, as discussed in Section II.2, follows a heavy tailed distribution with $\mu=5/2$ and $\Delta=1/2$ due to the nodal error in the trial wave function. In Fig. 13 we show $yx$ plots of the tails of the local Pulay force at sample sizes $M=10^{6}$, $10^{7}$, and $10^{8}$, along with plots of the corresponding “optimal” fits.

Despite having chosen a system known to exhibit a large nodal error, we find that the leading-order heavy tails are relatively weak, and that it takes sample sizes of $M\gtrsim 10^{7}$ to resolve the nonzero value of $y(0)$. As a result, the uncertainty on the standard estimator of $\langle A\rangle$ is likely to only exhibit nonconvergent behavior at large sample sizes.

We plot the convergence of the standard and tail-regression estimators of $\langle F_{\rm P}\rangle$ in Fig. 14.

As expected, the standard estimator $\bar{A}$ seems well behaved at small sample sizes, but at $M=10^{8}$ the standard error presents a substantial nonmonotonic jump. The uncertainty obtained with the tail-regression estimator remains smooth and monotonic, and is ultimately smaller than that of the standard estimator at $M=10^{8}$.

We find that the local zero-variance corrected Hellmann-Feynman force, $F_{\rm HFT}+F_{\rm ZV}$, and the local zero-variance corrected total force, $F_{\rm HFT}+F_{\rm P}+F_{\rm ZV}$, exhibit similarly weak leading-order tails, and the uncertainties in the standard and tail-regression estimators follow convergence patterns similar to those depicted in Fig. 14.

Without the zero-variance correction, the heavy tails affecting the distribution of the local Hellmann-Feynman force are very strong. As detailed in Section II.2, these tails are caused by the presence of all-electron nuclei and the left and right tails have equal leading-order coefficients. In Fig. 15 we plot the probability distribution of $F_{\rm HFT}$ for the off-equilibrium carbon dimer estimated from $10^{8}$ sample points and the corresponding $yx$ plots of the left and right tails of the distribution.

The value of $y(0)$ is very large relative to the rest of the function, and this is the only case among those we have considered in which the slope of the $yx$ plot is markedly negative at the origin. Indeed, it can be shown that the electron-nucleus cusp condition causes the $c_{1}$ coefficient in the asymptotic form of $F_{\rm HFT}$ to be approximately proportional to $c_{0}$ with a large, negative prefactor.

The standard and tail-regression estimators of the expectation value of the Hellmann-Feynman force are given in Table 5 and plotted in Fig. 16 as a function of sample size.

The uncertainty on the standard estimator is clearly nonmonotonic, while that in the tail-regression estimator is smooth and appears to present an $M^{-1/2}$ asymptotic decay. The uncertainty on the tail-regression estimator is up to $60$ times smaller than the nominal uncertainty on the standard estimator, of which a factor of $2.5$ is thanks to imposing the analytical constraint $c_{0}^{\rm L}=c_{0}^{\rm R}$.

The results obtained for the different force components at the largest considered sample size of $M=10^{8}$ are given in Table 6.

The nominal uncertainty on the standard estimator of $\langle F_{\rm HFT}+F_{\rm ZV}\rangle$ is an order of magnitude smaller than that in the constrained tail-regression estimator of $\langle F_{\rm HFT}\rangle$, and in this sense the tail-regression estimator is less effective than the zero-variance correction. This is however somewhat misleading since the uncertainty on the standard estimator remains formally ill-defined, while the tail-regression estimator is asymptotically normally distributed. In any case, the tail-regression estimator of $\langle F_{\rm HFT}+F_{\rm ZV}\rangle$ yields a $30\%$ lower uncertainty than the standard estimator, which is equivalent to a factor-of-two reduction in the number of sample points required to achieve a target uncertainty, showing that the combination of variance-reduction techniques with the tail-regression estimator is advantageous. Similarly, the nominal uncertainties on the Pulay force and on the zero-variance corrected total force are significantly reduced by replacing the ill-defined standard estimator with the tail-regression estimator at this sample size.

VII Application to diffusion quantum Monte Carlo data

At each post-equilibration step of a DMC calculation, an ensemble of walkers represents the mixed distribution $\Phi({\bf R})\Psi({\bf R})$, where $\Phi({\bf R})$ is the DMC wave function and $\Psi({\bf R})$ is the trial wave function. These walkers carry variable weights which in turn trigger death and branching events. The local values $\Psi^{-1}({\bf R})\hat{A}\Psi({\bf R})$ of an observable $\hat{A}$ are evaluated for each walker at each step, and the weighted average of the resulting sample yields the mixed estimator of the expectation value,

[TABLE]

Note that one in principle seeks the pure estimator,

[TABLE]

but the mixed estimator is simpler to obtain, it is equal to the pure estimator if $\hat{A}$ commutes with the Hamiltonian of the system, and for other observables there are ways of approximating pure estimators using mixed estimators Foulkes_QMC_2001 . We will restrict our discussion and tests to mixed estimators.

DMC samples differ from VMC samples in important ways. The formalism presented in section IV can be trivially altered to accommodate weights, simply by replacing the sample quantiles $q_{m}=\frac{m-1/2}{M}$ with

[TABLE]

where $p_{m}$ is the unnormalized weight of the $m$th sample point.

Walker branching events involve walkers being duplicated and its copies then evolving independently. This causes a complex pattern of serial correlation which cannot be eliminated entirely by leaving several steps between consecutive evaluations of the local values of observables, as we have done in our VMC calculations. While the presence of any form of serial correlation violates our assumption that samples consist of independent and identically distributed random variables, we expect this effect to be small and ignore it in our DMC tests.

The gradient of the DMC wave function $\Phi({\bf R})$ is in general discontinuous at the nodes Reynolds_dmc_1981 ; Moskowitz_lih_1982 . This alters the relative presence of walkers on either side of each nodal point, causing observables whose local values diverge at the nodes, such as the energy, to exhibit fully asymmetric heavy tails. The local Hellmann-Feynmann component of the force is not affected by this, since its singularities do not occur at the nodes of the trial wave function, so its DMC distribution remains symmetric to leading order as it is in VMC.

VII.1 The atomic force in the C2 molecule

We have performed a DMC simulation of the C2 molecule at the same off-equilibrium geometry and with the same wave function as described in Section VI.2, using a time step of $0.01$ a.u. Lee_strategies_2011 and a target population of $500$ walkers, and we have evaluated the local Hellmann-Feynmann and total forces for each walker every $5000$ steps.

The standard and tail-regression estimator of the Hellmann-Feynmann force are given in Table 7 and plotted in Fig. 17 as a function of sample size $M$. These results are very similar to their VMC counterparts; the nonmonotonic nominal standard error is again up to $60$ times the uncertainty in the tail-regression estimator, of which a factor of $2.5$ comes from imposing the constraint $c_{0}^{\rm L}=c_{0}^{\rm R}$.

The results we obtain for the total force at the largest considered sample size of $M=10^{8}$ are given in Table 8. In this case, the uncertainty in the tail-regression estimator of the total force is $10\%$ smaller than the standard error.

From these tests we conclude that the tail-regression estimator is directly applicable to DMC samples generated using relatively long decorrelation loops, with essentially the same benefits as we have found for VMC samples.

VIII Conclusions

We have introduced a conceptually simple estimator of expectation values for heavy tailed probability distributions whose power-law tail indices are known. Unlike the standard estimator, the tail-regression estimator is immune to the breakdown of the central limit theorem for distributions of leading-order tail exponent $2<\mu\leq 3$. Our regression framework is designed to yield asymptotically normally distributed results, as reflected in the observed asymptotic $M^{-1/2}$ decay with sample size $M$ of the uncertainty in all of our tests, and successfully exploits known analytical relations between leading order tail coefficients to improve the estimation. We have also demonstrated the estimation of the variance of distributions of leading-order tail exponent $3<\mu\leq 5$ whose uncertainty is ill-defined under standard estimation.

Our tests of the tail-regression estimator with variational and diffusion quantum Monte Carlo data identifies two use cases of particular practical relevance. While the standard estimator yields accurate expectation values of the energy at large enough sample sizes, standard confidence intervals on the VMC variance of the local energy are formally undefined. The tail-regression estimator is capable of delivering valid confidence intervals on the variance which are up to $45\%$ smaller than those associated with the nominal standard error in our tests. The tail-regression estimator also yields valid, convergent confidence intervals on the VMC and DMC atomic force, including the Hellmann-Feynman force in all-electron systems for which we obtain uncertainties up to $60$ times smaller than the nominal standard error. The combination of the “zero-variance” variance-reduction technique with the tail-regression estimator yields accurate confidence intervals on the atomic force.

Our present work shows that the principles underpinning the tail-regression estimator are robust, and systematic use of the technique for treating quantum Monte Carlo data would be desirable. However, further work could improve the applicability of our present formulation. We have used the bootstrap to enable the evaluation of meaningful confidence intervals on a range of functions, but this approach should be replaced with the use of a closed expression for the uncertainty on the tail-regression estimator in production calculations. In turn, this would allow the development of an “on-the-fly” reformulation of the method that would avoid the need to store all local values of the desired observables, unaveraged, for later analysis. Dropping the requirement that sample points be independent and serially uncorrelated would also be desirable in order to reduce the computational cost of the QMC calculations. With these refinements, the tail-regression estimator will ultimately represent a great advance in ensuring the statistical soundness of results obtained from quantum Monte Carlo and similar methods.

Acknowledgements.

The authors thank John Trail, Neil Drummond, and Richard Needs for useful discussions, and acknowledge the financial support of the Max-Planck-Gesellschaft, the Engineering and Physical Sciences Research Council of the United Kingdom under Grant No. EP/P034616/1, and the Royal Society. Supporting research data can be freely accessed at https://doi.org/10.17863/CAM.37836, in compliance with the applicable Open Data policies. Our implementation of the tail-regression estimator can also be found at treat_github and is distributed with the casino code casino_reference .

Appendix A Tail-index estimation methods

Tail-index estimation methods draw inference on the principal exponent $\mu$ of a power-law heavy tail of leading-order form $P_{A}(A)=c_{0}A^{-\mu}$ at $A\to\infty$. The Hill estimator Hill_TIE_1975 of the first-order tail index, $\mu-1$, is

[TABLE]

It can be shown that Eq. 48 is in fact equivalent to a logarithmic-scale least-squares fit to the tail of the distribution Beirlant_TIE_1996 . Substituting the leading-order form of $P_{A}(A)$ into Eq. 27 and taking logarithms yields

[TABLE]

which is a linear relationship between $\log A^{(m)}$ and $-\log q_{m}$ with slope $\frac{1}{\mu-1}$. Estimation of this slope by linear regression following Eq. 49 yields Eq. 48 for the fitted slope if the $m$th data point is weighted by

[TABLE]

and $c_{0}$ is set so that the fit passes through the $(M_{\rm R}+1)$th point. The optimal value of $M_{\rm R}$ for the Hill estimator can thus be found by optimizing a goodness-of-fit measure with respect to $M_{\rm R}$ Beirlant_TIE_1996 , such as the $\chi^{2}$ value of the fit. Regression methods for tail-index estimation are found to be particularly robust Baek_TIE_2010 .