text:
# Systematic Finite-Sampling Inaccuracy in Free Energy Differences and Other Nonlinear Quantities

## I Introduction

Because of the substantial recent interest in free energy difference $`\mathrm{\Delta }F`$ calculations (e.g., ), this report discusses the unavoidable error that arises from use of a finite amount of computer time. There is a tremendous range of applications for computational $`\mathrm{\Delta }F`$ estimates in physical, chemical, and biological systems. Examples include computations relating crystalline lattices , the behavior of magnetic models , and biomolecular binding events — of ligands to both DNA and proteins (e.g., ). Computations of $`\mathrm{\Delta }F`$, moreover, are formally equivalent to calculating the temperature dependence $`F(T)`$ . Most recently, it has been pointed out that $`\mathrm{\Delta }F`$ calculations are required to convert *experimental* data from nonequilibrium single-molecule pulling measurements to free energy vs. extension profiles ; see also . From a methodological standpoint, free energy computation protocols have been the subject of long and sustained interest

Free energy computations have long been recognized to suffer from a number of sources of error: (i) inaccuracy of the model (Hamiltonian) itself, (ii) incomplete conformational sampling, and (iii) finite sample size. In biomolecular systems, the issue of model accuracy (i) is indeed important, as typical all-atom force-fields are generally assumed to be capable of no more than 1 - 2 $`k_BT`$ accuracy in estimates of free energies — e.g., of ligand binding. That is, even with perfect sampling, computational estimates typically will not match experimental values. Second, like every simulation technique, free energy calculations are subject to errors based on (ii) incomplete conformational sampling. “*Incomplete*” sampling here refers to the lack of access to important parts of phase or conformational space: that is, the distribution of samples of size $`N`$ generated by the simulation will *not* match the true (representative) distribution of size-$`N`$ samples which would be drawn at random from a very long, perfectly-sampled simulation. Incomplete conformational sampling introduces bias into even simple computations attempting to estimate linear quantities, such as the mean of some coordinate or function. Conformational sampling errors in free energy calculations have long been recognized as “hysteresis” or “Hamiltonian-lag,” and a number of workers have made important contributions toward understanding and overcoming these errors — e.g., .

We focus here on the third type of error (iii) that due *solely* to the necessarily finite samples collected in a simulation. Such finite-sampling bias occurs in every non-linear calculation, as detailed below, and should be clearly distinguished from the independent error due to (ii) inadequate conformational sampling. Specifically, finite-sampling error occurs even when conformational sampling is perfect — i.e., when representative samples of finite size are generated by the simulation. Finite-sampling errors in computational estimates of $`\mathrm{\Delta }F`$ were first recognized by Wood and coworkers  and later discussed by others ; see also . The in-depth work of Lu and Kofke presents an entropy-based description of finite-sampling errors, which emphasizes the critical asymmetry between generalized “insertion” and “deletion” calculations .

Figure 1 illustrates the phenomenon of finite-sampling error for a mathematical model and for a biological system , emphasizing the universality of finite-sampling errors. Because these inaccuracies can be many times $`k_BT`$ (see Fig. 1 and Ref. ) — especially in the important context of biomolecular calculations where large system sizes limit the quantity of data available for analysis — there is a strong motivation to understand and overcome these errors. Ferrenberg, Landau and Binder showed that finite-sampling errors accompanying susceptibility computations can be understood on the basis of elementary statistical principles ; however, the errors in *non-linear averages* like $`\mathrm{\Delta }F`$ apparently had remained without an explicit theoretical basis until recently . In a recent monograph, in fact, Landau and Binder note that finite-sampling errors are “generally given inadequate attention” .

This report both provides fuller details of the theory outlined in , and also presents new results. Our report includes (i) a detailed proof that the expected value of a finite-data $`\mathrm{\Delta }F`$ estimate ($`\mathrm{\Delta }F_n`$) bounds the true free energy — *independent of the distribution of underlying work values*; (ii) full derivations of the asymptotic expressions for $`\mathrm{\Delta }F_n`$ as $`n\mathrm{}`$ for arbitrary — including long-tailed — distributions of the work ($`W`$) data used to estimate $`\mathrm{\Delta }F`$; (iii) analogous derivations for the root-mean-square and related “geometric” non-linear averages; (iv) derivation and numerical demonstration of the universal asymptotic relation between $`\mathrm{\Delta }F_n`$ and its fluctuation. As in our brief report , the present discussion makes use of mathematical results regarding the convergence — to “stable” limiting distributions , also known as Lévy processes (e.g., ) — of the distributions of sums of variables.

In outline, the paper now proceeds to Sec. II where formal groundwork for the discussion is laid. Section III rigorously proves the true free energy $`\mathrm{\Delta }F`$ is bounded by $`\mathrm{\Delta }F_n`$, the expected value of a free-energy estimate based on a finite quantity of data; analogous bounds apply for *arbitrary* non-linear averages. Derivations of the asymptotic series for $`\mathrm{\Delta }F_n`$ are given in sections IV and V, while Section VI derives the universal relation between $`\mathrm{\Delta }F_n`$ and its fluctuation. We conclude with a summary and discussion of the results in Section VII.

## II Free-Energy Estimates from Finite Samples

Since the work of Kirkwood , it has been appreciated that the free energy difference, $`\mathrm{\Delta }F\mathrm{\Delta }F_{01}`$, of switching from a Hamiltonian $`_0`$ to $`_1`$ is given by a non-linear average,

$$\mathrm{\Delta }F=k_BT\mathrm{log}\left[\mathrm{exp}(W_{01}/k_BT)_0\right],$$
(1)
where $`k_BT`$ is the thermal unit of energy at temperature $`T`$ and $`W_{01}`$ is the work required to switch the system from $`_0`$ to $`_1`$; see below. The angled brackets indicate an average over switches starting from configurations drawn from the equilibrium distribution governed by $`_0`$.

While non-equilibrium approaches to free energy calculations have been a major motivation for this work, we should point out that our analysis applies equally to “staged” calculations, in which the free energy is calculated as a sum of increments. In particular, if one writes the free energy as a sum of incremental parts,

$$\mathrm{\Delta }F_{01}=\mathrm{\Delta }F_{0\lambda _1}+\mathrm{\Delta }F_{\lambda _1\lambda _2}+\mathrm{}+\mathrm{\Delta }F_{\lambda _k1},$$
(2)
then each increment $`\mathrm{\Delta }F_{\lambda _i\lambda _j}`$ is still defined by a non-linear average analogous to (1) and thus will suffer from finite-sampling error.

The work $`W_{01}`$ required for the average (1) can be defined in a straightforward manner. In instantaneous switching the work is defined by $`W_{01}=_1(\text{x})_0(\text{x})`$ for a start (and end) configuration x. However, gradual switches requiring a “trajectory”-based work definition may also be used, as was demonstrated by Jarzynski . In this latter case, one requires a Hamiltonian which interpolates between $`_0`$ and $`_1`$; a common choice is

$$(\lambda ;\text{x})_0(\text{x})+\lambda \left[_1(\text{x})_0(\text{x})\right],$$
(3)
where $`\lambda `$ is a switching parameter that varies from zero to one. The work performed in switching gradually from $`_0`$ to $`_1`$ along a trajectory $`(\lambda (t);\text{x}(t))`$ is given by

$$W_{01}=\underset{i}{}\left[(\lambda _i;\text{x}_{i1})(\lambda _{i1};\text{x}_{i1})\right],$$
(4)
where the subscripted configuration $`\text{x}_{i1}`$ is the (unique) *final* configuration for which $`\lambda =\lambda _{i1}`$ — i.e., the last configuration before $`\lambda `$ is incremented to $`\lambda _i`$. In other words, the work is computed as the sum of those energy increments resulting only from changes in $`\lambda `$.

Whenever a convex, nonlinear average such as (1) is estimated computationally, that result will *always* be systematically biased  because one has only a finite amount of data — say, $`N`$ work values. The bias results from incomplete sampling of the smallest (or most negative) $`W_{01}`$ values : these values dominate the average (1) and cannot be sampled perfectly for finite $`N`$, regardless of the $`W_{01}`$ distribution. This is true even for a rectangular distribution; the sole exception is the trivial $`\delta `$ function, single-point probability density. Because of the undersampling of small work values, a running estimate of $`\mathrm{\Delta }F`$ will typically decline as data is gathered, as one sees in the “staircase” plots of Fig. 1. Such considerations led Wood *et al.*  to consider the block-averaged $`n`$-data-point estimate of the free energy based on $`N=mn`$ total work values $`\{W^{(k)}\}`$, namely,

$$\mathrm{\Delta }F_n=\underset{m\mathrm{}}{lim}\frac{1}{m}\underset{j=1}{\overset{m}{}}k_BT\mathrm{log}\left[\frac{1}{n}\underset{k=(j1)n+1}{\overset{jn}{}}\mathrm{exp}(W^{(k)}/k_BT)\right].$$
(5)
It represents the expected value (mean) of a free energy estimate from $`n`$ data points — that is, of

$$_n=k_BT\mathrm{log}\left[\left(e^{W_1/k_BT}+\mathrm{}+e^{W_n/k_BT}\right)/n\right],$$
(6)
where $`m`$ estimates have been made. Wood *et al.* estimated the lowest order correction to $`\mathrm{\Delta }F\mathrm{\Delta }F_{\mathrm{}}`$ as $`\sigma _w^2/2nk_BT`$, where $`\sigma _w^2`$ is the variance in the distribution of work values, $`W`$ . Ferrenberg, Landau and Binder discussed analogous issues for the magnetic susceptibility .

The derivations below employ continuum expressions simplified by the definitions

$$wW/k_BT,f\mathrm{\Delta }F/k_BT,f_n\mathrm{\Delta }F_n/k_BT.$$
(7)
In terms of the probability density $`\rho _w`$ of work values, which is normalized by $`dw\rho _w(w)=1`$, the free energy is given by the continuum analog of (1),

$$f=\mathrm{\Delta }F/k_BT=\mathrm{log}\left[dw\rho _w(w)e^w\right].$$
(8)
The form (8) also occurs in equilibrium calculations , and forms the basis for the analysis of Lu and Kofke . Finally, the finite-data average free energy, following (5) must apply the logarithm “before” the average of the $`n`$ Boltzmann factors, and one has

$`f_n`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}\left[\mathrm{d}w_i\rho _w(w_i)\right]\mathrm{log}\left[\frac{1}{n}\underset{i=1}{\overset{n}{}}e^{w_i}\right]}.`$ (9)
## III Generalized Jensen’s Inequalities: Block-Averaged Estimates as Rigorous Bounds

Jarzynski  and subsequently the present authors  observed that the free energy is bounded according to

$$\mathrm{\Delta }F\mathrm{\Delta }F_n,\text{any }n.$$
(10)
Here we prove this inequality and a generalization originally stated in , namely,

$$\mathrm{\Delta }F_{n+1}\mathrm{\Delta }F_n,\text{any }n.$$
(11)
Indeed the proof given below applies to a broad class of nonlinear averages and functions. As noted in , the result (10) substantially extends the previous bound $`\mathrm{\Delta }FW\mathrm{\Delta }F_1`$ . Further, the bounds apply *for an arbitrary distribution of work values* — that is, whether the probability density of $`W`$ is multimodal, unimodal, or even rectangular. We note, finally, that reference  by the present authors failed to acknowledge the original statement of the bound (10) by Jarzynski in .

In fact, our proof will show that (10) and (11) are special cases of a more general inequality that depends solely on the convexity and monotonicity of the function used to form a nonlinear average: in the case of $`\mathrm{\Delta }F`$ the function is the exponential — see (1); the root-mean-square is another example, when the function is $`g(x)=x^2`$. In the remainder of this section, we use the mathematical convention that upper-case letters (e.g., $`X`$) indicate random variables whose particular values are specified by lower-case letters (e.g., $`x`$).

The new bounds are generalizations of Jensen’s inequality (see ), a fundamental property of convex functions with a host of applications including in information theory . Jensen’s inequality relates the expected value of a convex function $`g`$ of a random variable to the same function of the expected value of its argument according to

$$g(X)g(X),$$
(12)
where the expectation value is defined in the usual way for an arbitrary function $`A`$ as

$$A=A(X)=dx\rho (x)A(x),$$
(13)
and $`\rho `$ is the probability density function. By applying $`g^1`$ to (12), the inequality can be re-stated in terms of non-linear and linear averages, respectively,

$$X^gg^1\left(g(X)\right)X,$$
(14)
with the additional constraint that $`g`$ be increasing (so that $`g^1`$ is unique). Note that the inequality (14) can easily be generalized by applying the inverse of a different increasing function (say, $`h^1`$) to (12).

We now state and prove the new “generalized Jensen’s inequalities.”
*Theorem:*
Consider estimates for the non-linear average $`X^g`$ based on samples of size $`n`$, $`\{x_1,x_2,\mathrm{},x_n\}`$, the expectation of which may be written as

$$X_n^g=dx_1\rho (x_1)\mathrm{}dx_n\rho (x_n)g^1\text{(}\left[g(x_1)+\mathrm{}g(x_n)\right]/n\text{)}.$$
(15)
Note that $`X_1^g=X`$ and $`X_{\mathrm{}}^g=X^g`$. Then the new inequalities, generalizing (14), are

$$X_n^gX_{n1}^g.$$
(16)
Strict inequality obtains whenever the random variable $`X`$ is not limited to a single value (i.e., whenever the probability density $`\rho `$ is not a single Dirac delta function). The direction of the inequality is reversed for *decreasing* convex functions, for instance yielding (10) for $`g(x)=\mathrm{exp}(x/k_BT)`$.

*Proof:*
Note first that $`X_n^g`$ is defined in (15) as the non-linear average based on the “weighted set” $`S_n`$ of all possible $`n`$-samples $`\{x_1,\mathrm{},x_n\}`$. The weight of each $`n`$-sample is of course its probability density $`\rho _n(\{x_i\})=_{i=1}^n\rho (x_i)`$. We will require an explicit construction of the set $`S_{n1}`$ from $`S_n`$, which fortunately is straightforward: for every $`n`$-sample with weight $`\rho _n`$ in $`S_n`$, if one assigns equal weights $`\rho _n/n`$ to each of the $`n`$ available ($`n1`$)-samples given by deletion of a single element — namely, $`\{x_2,x_3,\mathrm{},x_n\}`$, $`\{x_1,x_3,x_4,\mathrm{},x_n\}`$, and so on — one arrives at $`S_{n1}`$. The correctness of this construction follows from the factorizability of the density $`\rho _n`$, and may be seen by considering the density of a particular ($`n1`$)-sample, $`\{\widehat{x}\}=\{\widehat{x}_1,\mathrm{},\widehat{x}_{n1}\}`$, which can be constructed from $`n`$ different deletions:

$`\rho _{n1}(\{\widehat{x}\})`$ $`=`$ $`{\displaystyle \frac{1}{n}}[{\displaystyle }\mathrm{d}x\rho _n(x,\widehat{x}_1,\mathrm{},\widehat{x}_{n1})+{\displaystyle }\mathrm{d}x\rho _n(\widehat{x}_1,x,\widehat{x}_2,\mathrm{},\widehat{x}_{n1})`$ (18)
$`+\mathrm{}+{\displaystyle }\mathrm{d}x\rho _n(\widehat{x}_1,\mathrm{},\widehat{x}_{n1},x)]`$
$`=`$ $`{\displaystyle \frac{1}{n}}\left[n{\displaystyle dx\rho (x)\underset{i=1}{\overset{n1}{}}\rho (\widehat{x}_i)}\right]={\displaystyle \underset{i=1}{\overset{n1}{}}}\rho (\widehat{x}_i).`$ (19)
Because of this construction of $`S_{n1}`$ from $`S_n`$, it is sufficient to show that the *single-sample* non-linear average of an arbitrary $`n`$-sample, namely,

$$u_n(\{x_i\})=g^1\left(\frac{1}{n}\underset{i=1}{\overset{n}{}}g(x_i)\right)$$
(20)
exceeds the average $`u_{n1}`$ based on the $`n`$ available ($`n1`$)-samples constructed from deletions, as above. Note that

$$u_n(\{x_i\})_n=X_n^g,$$
(21)
which follows immediately from (15).

To complete the proof, observe that the single-sample non-linear average (20) can be re-written in terms of smaller samples:

$$u_n(\{x_i\})=g^1\left(\frac{1}{n}\underset{j=1}{\overset{n}{}}\frac{1}{n1}\underset{ij}{\overset{n}{}}g(x_i)\right).$$
(22)
This identity may be illustrated by considering $`g(x_1)`$ which occurs $`n1`$ times (whenever $`j1`$), and hence is properly weighted as in (20). The expression may be further re-written if we denote by $`\{x_i\}_{[j]}`$ the original $`n`$-sample with the $`j`$th element deleted. To each of these smaller samples, there corresponds a single-sample, non-linear average $`u_{n1}(\{x_i\}_{[j]})`$. Applying $`g`$ to both sides of (20) and substituting the result for $`n1`$ into the right-hand-side of (22), we then have

$$u_n(\{x_i\})=g^1\left(\frac{1}{n}\underset{j=1}{\overset{n}{}}g\left[u_{n1}\text{(}\{x_i\}_{[j]}\text{)}\right]\right).$$
(23)
If we now consider $`U_{n1}u_{n1}`$ to be a random variable with a discrete, $`n`$-point distribution, we can apply the original non-linear-average inequality (14) to the right-hand side of (23), and obtain the desired result

$$u_n(\{x_i\})u_{n1}_{[j]},$$
(24)
where the average $`\mathrm{}_{[j]}`$ is performed over the discrete distribution comprised of all $`u_{n1}`$ values obtained from applying (20) to the $`n`$ sets $`\{x_i\}_{[j]}`$. This completes the proof because when the left-hand-side is averaged as in (21), the construction of $`S_{n1}`$ from $`S_n`$ guarantees that the average over all $`n`$-samples on the right-hand side of (24) results in $`X_{n1}^g`$ and hence (16).

The result applies to *any* probability density $`\rho `$ because no assumptions were made regarding the distribution. Figure 1 illustrates the monotonicity of finite-data free energy estimates from two completely unrelated systems.

## IV Asymptotic Behavior: Finite Moments Case

### A Formal Development of the Expansion

It is possible to generate a formal expansion for the finite-data estimate $`f_n`$ in terms of $`n^1`$ for an arbitrary distribution of work values $`\rho _w`$. In this section we consider the case where the second and some higher moments of the $`z=e^w`$ distribution are finite. Motivated by the central and related limit theorems  for the sum of the $`e^w`$ variables, we introduce a change of variables which will permit the development of a $`1/n`$ expansion for $`f_n`$. In particular, we define

$$y=(e^{w_1}+\mathrm{}+e^{w_n}ne^f)/b_1n^{1/\alpha },$$
(25)
where $`b_1`$ is a constant and $`\alpha 2`$ is an exponent characterizing the distribution of the variable $`e^w`$. The requirement that $`\mathrm{\Delta }F`$ be finite in (8) further implies $`\alpha >1`$. The finite-data free energy difference can now be written

$$f_n=_{cn^a}^{\mathrm{}}dy\rho _n(y)\mathrm{log}\left(e^f+\frac{b_1}{n^a}y\right),$$
(26)
where $`c=\mathrm{exp}(f)/b_1`$, $`a(\alpha 1)/\alpha <1/2`$, and $`\rho _n`$ is the probability density of the variable $`y`$ normalized appropriately via

$$_{cn^a}^{\mathrm{}}dy\rho _n(y)=1.$$
(27)
Note that $`a`$ is always positive because $`\alpha >1`$.

The expansion of $`f_n`$ proceeds by first noting that the sum of *any* set of independent random variables, suitably normalized as in (25), has a distribution which may be expressed as a stable (Lévy) distribution function multiplied by a large-$`n`$ asymptotic expansion . In the case of a Gaussian limiting distribution (i.e., $`\alpha =2`$ or the central limit theorem), assume the variable $`z=e^w`$ possesses finite “Boltzmann moments” — a mean $`\widehat{\mu }=e^f`$, variance $`\widehat{\sigma }^2=(z\widehat{\mu })^2`$, and higher central moments $`\widehat{\mu }_p=(z\widehat{\mu })^p`$. The normalizing constant in (25) is then $`b_1=\widehat{\sigma }`$. The Boltzmann moments of course differ from the moments of the distribution of $`w`$.

The so-called Edgeworth corrections to the central limit theorem indicate that the variable $`y=(_{i=1}^ne^{w_i}n\widehat{\mu })/\sqrt{n}\widehat{\sigma }`$ \[cf. (25)\] is distributed according to

$$\rho _n(y)=\rho _G(y;1)\left[1+\nu _1(y)/\sqrt{n}+\nu _2(y)/n+\mathrm{}\right],$$
(28)
for large $`n`$, where the remaining terms are higher integer powers of $`1/\sqrt{n}`$ and the Gaussian density is

$$\rho _G(y;\sigma )=\mathrm{exp}(y^2/2\sigma ^2)/\sqrt{2\pi }\sigma .$$
(29)
The functions $`\nu _i`$, which are defined based upon the Hermite polynomials , depend on the original distribution of $`e^w`$. In terms of the cumulants $`\widehat{\kappa }_i`$ (see, e.g., ) of the distribution of $`z=e^w`$ and the Hermite polynomials defined via

$$\frac{\mathrm{d}^k}{\mathrm{d}x^k}\rho _G(x;1)=(1)^kH_k(x)\rho _G(x;1),$$
(30)
the lowest-order Edgeworth functions are

$`\nu _1(y)`$ $`=`$ $`(\widehat{\kappa }_3/6\widehat{\sigma }^3)H_3(y)=(\widehat{\mu }_3/6\widehat{\sigma }^3)\left(y^33y\right)`$ (31)
$`\nu _2(y)`$ $`=`$ $`(\widehat{\kappa }_4/24\widehat{\sigma }^4)H_4(y)+(\widehat{\kappa }_3^2/72\widehat{\sigma }^6)H_6(y)`$ (32)
$`\nu _3(y)`$ $`=`$ $`(\widehat{\kappa }_5/120\widehat{\sigma }^5)H_5(y)+(\widehat{\kappa }_3\widehat{\kappa }_4/144\widehat{\sigma }^7)H_7(y)+(\widehat{\kappa }_3^3/1296\widehat{\sigma }^9)H_9(y).`$ (33)
The $`\nu _i`$ functions are odd or even according to whether $`i`$ is odd or even, in this $`\alpha =2`$ case.

Before the expansion for $`f_n`$ can be developed, the integral (26) must be considered carefully by dividing it into three parts:

$`f_n`$ $`=`$ $`{\displaystyle _{cn^a}^{\mathrm{}}}dy\rho _n(y)\mathrm{log}\left(e^f+{\displaystyle \frac{b_1}{n^a}}y\right)`$ (34)
$`=`$ $`{\displaystyle _{cn^a}^{\mathrm{}}}dy\rho _n(y)\mathrm{log}\left(e^f\right)`$ (36)
$`+{\displaystyle _{cn^a}^{+cn^a}}dy\rho _n(y)\mathrm{log}\left(1+y/cn^a\right)+{\displaystyle _{+cn^a}^{\mathrm{}}}dy\rho _n(y)\mathrm{log}\left(1+y/cn^a\right)`$
$``$ $`f+I(cn^a,cn^a)+I(cn^a,\mathrm{}),`$ (37)
where the first integral in (36) has been evaluated exactly using the normalization of $`\rho _n`$ (27) and $`I`$ represents the latter integrals of (36). One can now proceed by using an expansion for the logarithm in $`I(cn^a,cn^a)`$ and by bounding terms in $`I(cn^a,\mathrm{})`$.

It is possible to demonstrate rigorously that the second integral in (36), $`I(cn^a,\mathrm{})`$, does not materially contribute to $`f_nf`$ for large $`n`$. Although, the logarithm cannot be expanded in a power series for $`y>cn^a`$, the integral can be bounded by expressing the log as the integral of its derivative:

$`I(cn^a,\mathrm{})`$ $`=`$ $`{\displaystyle _{+cn^a}^{\mathrm{}}}dy\rho _n(y){\displaystyle _1^{1+y/cn^a}}dxx^1`$ (38)
$``$ $`{\displaystyle _{+cn^a}^{\mathrm{}}}dy\rho _n(y){\displaystyle _1^{1+y/cn^a}}dxx^{1+ϵ}={\displaystyle \frac{1}{ϵ}}{\displaystyle _{+cn^a}^{\mathrm{}}}dy\rho _n(y)\left[\left(1+{\displaystyle \frac{y}{cn^a}}\right)^ϵ1\right],`$ (39)
with $`0<ϵ1`$. To extract the leading behavior of this bound, one can use the expansion of $`\rho _n`$ (28) and set $`ϵ=1`$. Noting that $`a=1/2`$, one obtains

$$I(cn^a,\mathrm{})\frac{1}{c\sqrt{n}}_{+cn^a}^{\mathrm{}}dy\rho _G(y;1)\left[1+\nu _1(y)/\sqrt{n}+\nu _2(y)/n+\mathrm{}\right]y.$$
(40)
Using the asymptotic properties of the error function , one can show that the strongest $`n`$ dependence of $`I(cn^a,\mathrm{})`$ is no stronger than

$$n\mathrm{exp}(c^2n/2).$$
(41)
The leading behavior of $`f_nf`$ is thus expected to result from the first non-trivial integral in (36), $`I(cn^a,cn^a)`$. Noting again that $`a=1/2`$ in this case, we may write

$`I(c\sqrt{n},c\sqrt{n})`$ $`=`$ $`{\displaystyle _{c\sqrt{n}}^{+c\sqrt{n}}}dy\rho _G(y;1)\left[1+\nu _1(y)/\sqrt{n}+\nu _2(y)/n+\nu _3(y)/n^{3/2}+\mathrm{}\right]`$ (43)
$`\times \left[y/c\sqrt{n}(y/c\sqrt{n})^2/2+(y/c\sqrt{n})^3/3\mathrm{}\right].`$
What are the leading terms? There are no terms proportional to $`n^{1/2}`$ raised to any odd power because of symmetry considerations: the $`\nu _i`$ functions are even for even $`i`$. The leading terms are thus *integer powers* of $`n^1`$, and the expansion of the finite-data free-energy difference is of the form,

$$f_n=f+\phi _1/n+\phi _2/n^2+\mathrm{},$$
(44)
where the $`\phi _i`$ are constants which depend on the distribution of $`z=e^w`$.

The explicit correction terms to $`f_nf`$ may now be obtained. First note that asympototic analysis of the integrals appearing in (43) in terms of the error function  indicates that the limits of integration may be extended to $`(\mathrm{},+\mathrm{})`$ with errors proportional to $`\mathrm{exp}(c^2n/2)`$. Straightforward integration then yields the coefficients of the expansion (44), namely,

$`\phi _1`$ $`=`$ $`\widehat{\sigma }^2/2\widehat{\mu }^2,`$ (45)
$`\phi _2`$ $`=`$ $`(4\widehat{\mu }\widehat{\mu }_39\widehat{\sigma }^4)/12\widehat{\mu }^4.`$ (46)
### B Coefficients for the Gaussian case

When the distribution of work values is Gaussian, $`\rho _w(W)=\rho _G(W,\sigma _w)`$, the Boltzmann moments and, hence the $`\phi `$ coefficients of (44), may be computed analytically. Note that one *cannot* assume that $`z=e^w`$ obeys a Gaussian distribution because $`z`$ is always non-negative. The moments follow from straightforward integration, which yields

$$z^p=dW\rho _w(W)e^{pW/k_BT}=\mathrm{exp}\left[p^2\sigma _w^2/\mathrm{\hspace{0.17em}2}(k_BT)^2\right].$$
(47)
The $`f_n`$ expansion coefficients then follow trivially from substitution into (45) and (46). Setting $`s=\sigma _w/k_BT`$, one finds for the first two coefficents

$`\phi _1`$ $`=`$ $`\left(e^{s^2}1\right)/2,`$ (48)
$`\phi _2`$ $`=`$ $`\left(4e^{3s^2}+9e^{2s^2}6e^{s^2}+1\right)/12.`$ (49)
To compare this with the finding of Wood *et al.* for $`f_nf`$, one can expand (48) for small $`\sigma _w`$. One finds $`\phi _1\sigma _w^2/2k_BT`$, which is precisely the first-order prediction of Wood *et al.* .

This analytic calculation explicitly indicates the practical shortcomings of the expansion (44). Although the leading term in $`f_nf`$ is linear in $`1/n`$, the leading coefficients are exponential in the *square* of the distribution’s width. The asymptotic expression (44) thus represents a viable approximation only for a very small window about $`1/n=0`$ when $`s1`$; see Fig. 2. When asymmetry is added to a Gaussian distribution via the first Edgeworth correction (see (28) and, e.g., ), one finds that the exponential dependence of the $`\phi _i`$ on $`\sigma _w`$ is only corrected linearly by the now non-zero third moment of the $`W`$ distribution.

### C Expansions for the root-mean-square and similar averages

The root-mean-square (or standard deviation) is perhaps the best known example of a non-linear average. The full analysis carried out above carries over quite directly, and indeed applies to any non-linear average. We will now briefly consider general “root-mean-powers” (“power means”).

To be specific, consider the non-linear average resulting from a general power $`q=2,4,\mathrm{}`$, denoted

$$R^{(q)}x^q^{1/q},$$
(50)
where $`x`$ is a variable distributed according to the (arbitrary) probability density $`\rho _x`$.

In direct analogy with (9) one can define the finite-data average for $`R^{(q)}`$ as

$$R_n^{(q)}=\underset{i=1}{\overset{n}{}}\left[\mathrm{d}x_i\rho _x(x_i)\right]\left[\frac{1}{n}\underset{i=1}{\overset{n}{}}x_i^q\right]^{1/q}.$$
(51)
The asymptotic expansion follows from the same procedure as above. One finds that the expansion

$$R_n^{(q)}=R^{(q)}+\phi _1^{(q)}/n+\phi _2^{(q)}/n^2+\mathrm{},$$
(52)
has coefficients

$`\phi _1^{(q)}`$ $`=`$ $`(1/2q)(1q^1)R^{(q)}\stackrel{~}{\sigma }^2/\stackrel{~}{\mu }^2`$ (53)
$`\phi _2^{(q)}`$ $`=`$ $`\left[q^1(1q^1)(2q^1)R^{(q)}/\stackrel{~}{\mu }^4\right]\left[\stackrel{~}{\mu }\stackrel{~}{\mu }_3/6(3q^1)\stackrel{~}{\sigma }^4/8\right]`$ (54)
where $`\stackrel{~}{\mu },\stackrel{~}{\sigma }^2,\text{ and }\stackrel{~}{\mu }_3`$ denote the mean, variance, and third central moment — respectively — of the distribution of the variable $`x^q`$.

## V Asymptotic Behavior: Divergent Moments Case

When the distribution $`\rho _z`$ of the variable $`e^wz`$ in (25) possesses a long-tail, the limiting distribution is not a Gaussian and the results (28) and (44) no longer hold. In particular, if one of the tails of $`\rho _z(z)`$ decays as $`z^{(1+\alpha )}`$ with $`\alpha <2`$ (implying an infinite Boltzmann variance, $`\widehat{\sigma }^2`$), then the distribution of the variable $`y`$ in (25) approaches a non-Gaussian “stable” (Lévy) law for large $`n`$ . Note that such power-law behavior in $`z`$ corresponds to *simple exponential decay in the work distribution.*

A long-tailed $`z`$ distribution $`\rho _z\rho _1`$ also alters the *form* of the asymptotic expansion of the distribution of the sum-variable (25) and, hence, the expansion of $`f_n`$ — which no longer includes solely integer powers of $`n^1`$, as in (44). Instead of (28), the $`y`$ distribution now takes the more complicated form

$$\rho _n(y)=\rho _\alpha (y)\left[1+\stackrel{}{}\nu _{uv}(y)/n^{\theta (u,v)}\right],$$
(55)
where $`\rho _\alpha `$ is the appropriate stable probability density with exponent $`\alpha `$ . The functions $`\nu _{uv}`$, which are not available analytically, depend on the original distribution of $`e^w`$ and partial derivatives of the stable distribution. The exponents are given by $`\theta (u,v)=(u+\alpha v)/\alpha >0`$, and the summation $`^{}`$ includes $`u0`$ and $`vu/2`$, where $`x`$ denotes the integer part of $`x`$.

To analyze the asymptotic behavior of $`f_n`$ in this case, the starting point is again equations (25) - (27), which are fully general. It is useful to rewrite (26) by scaling the logarithm’s argument by the constant $`e^f`$ and by subtracting zero in the form of the mean of $`y`$; one obtains

$$ff_n=_{cn^a}^{\mathrm{}}dy\rho _n(y)\left[\mathrm{log}\left(1+\frac{y}{cn^a}\right)\frac{y}{cn^a}\right]\widehat{I}(cn^a,\mathrm{}).$$
(56)
One can now divide up the domain of integration in (56) into sub-parts appropriate for expansions of the logarithm of $`\rho _n`$, in analogy with (36). Because no explicit forms for stable distributions are known in the range $`1<\alpha <2`$ , we will require separate expansions of $`\rho _n`$ for $`|y|1`$ and $`y\pm \mathrm{}`$ to obtain appropriate convergent behavior. The required breakdown of the integral is therefore

$$ff_n=\widehat{I}(cn^a,1)+\widehat{I}(1,1)+\widehat{I}(1,\mathrm{}).$$
(57)
Each of the integrals in (57) requires a slightly different procedure. The first, $`\widehat{I}(cn^a,1)`$, requires an expansion of the logarithm along with the “short-tail” $`y\mathrm{}`$ expansion of $`\rho _n\rho _\alpha `$ (see below). The second integral, $`\widehat{I}(1,1)`$, uses simple convergent series expansions of both the logarithm and $`\rho _\alpha `$. Finally, $`\widehat{I}(1,\mathrm{})`$ requires primarily the “long-tail” $`y\mathrm{}`$ expansion of $`\rho _\alpha `$; the series expansion of the logarithm is also used to show that extending the lower limit of integration to zero accrues a non-leading correction.

Because we will extract only the leading term of $`f_nf`$, it is sufficient to use only the leading contribution to $`\rho _n`$; that is, considering (55) we may use the asymptotically valid ($`n\mathrm{}`$) approximation $`\rho _n\rho _\alpha `$. (The leading behavior for $`f_n`$ in the finite-moments case arises, similarly, from $`\rho _n\rho _G`$.) The required series expansions for $`\rho _\alpha `$ in the case of positive summands $`z=e^w`$ are

$`\rho _\alpha (y;\xi )`$ $`={\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}C_k^0|y|^{k1},`$ $`|y|>0`$ (58)
$`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}C_k^{\mathrm{}}y^{(k\alpha +1)},`$ $`y\mathrm{}`$ (59)
where the “$``$” sign denotes an asympotic expansion, and the coefficients — which depend on the sign of $`y`$ — are given by

$`C_k^0(\xi )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}(1)^{k1}{\displaystyle \frac{\mathrm{\Gamma }(1+k/\alpha )}{k!}}\mathrm{sin}(k\pi \xi /\alpha ),`$ (60)
$`C_k^{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}(1)^{k1}{\displaystyle \frac{\mathrm{\Gamma }(k\alpha +1)}{k!}}\mathrm{sin}(k\pi \xi ^+),`$ (61)
with $`\xi ^+\xi (y>0)=\alpha 1`$ and $`\xi ^{}\xi (y<0)=1`$. Note that $`C_k^0(\xi ^+)=(1)^{k1}C_k^0(\xi ^{})`$, and in particular, $`C_1^0(\xi ^+)=C_1^0(\xi ^{})C_1^0`$. Because the summands considered here are strictly positive, the left tail of the distribution does not exhibit power-law behavior; rather, it may be termed “short” or “light” and, asymptotically, is given by

$$\rho _\alpha (y\mathrm{})\frac{1}{\sqrt{2\pi \alpha (\alpha 1)}}\left|\frac{y}{\alpha }\right|^{\frac{1\alpha /2}{\alpha 1}}\mathrm{exp}\left\{(\alpha 1)\left|\frac{y}{\alpha }\right|^{\frac{\alpha }{\alpha 1}}\right\}.$$
(62)
We can now consider the terms in (57) using (58) - (62). For the sake of brevity we quote only the leading terms, which result from straightforward integrations (after discarding non-leading terms and corrections):

$`\widehat{I}(cn^a,1)`$ $``$ $`{\displaystyle \frac{\alpha ^2\mathrm{\Gamma }(a+\frac{3}{2},\alpha ^a(\alpha 1))}{2c^2(\alpha 1)^{a+\frac{1}{2}}\sqrt{2\pi \alpha (\alpha 1)}}}n^{2a}`$ (63)
$`\widehat{I}(1,1)`$ $``$ $`\left(C_1^0/3c^2\right)n^{2a}`$ (64)
$`\widehat{I}(1,\mathrm{})`$ $``$ $`\left(C_1^{\mathrm{}}I_\alpha ^{}/c^\alpha \right)n^{1\alpha },`$ (65)
where $`2a=2(\alpha 1)/\alpha `$, $`\mathrm{\Gamma }(,)`$ is the incomplete gamma function , and

$$I_\alpha ^{}=_0^{\mathrm{}}\frac{\mathrm{d}x}{x^{1+\alpha }}\left[x\mathrm{log}(1+x)\right]<\mathrm{}.$$
(66)
By comparing powers of $`n`$ in (63) - (65) one sees that the leading behavior of the finite-date free energy estimate, not surprisingly, results from the “heavy” power-law tail ($`y\mathrm{}`$). Thus, using (65), one has

$$f_nf\phi _{\alpha 1}/n^{(\alpha 1)},$$
(67)
with $`\phi _{\alpha 1}=\left(C_1^{\mathrm{}}I_\alpha ^{}/c^\alpha \right)>0`$. Note that $`\phi _{\alpha 1}`$ depends on $`\alpha `$ and also on the original probability density $`\rho _z`$ through $`c=e^f/b_1`$. Furthermore, one should not expect (67) to be a useful estimate for $`f_nf`$: the next leading exponent, $`2(\alpha 1)/\alpha `$ is very close to $`\alpha 1`$ for $`\alpha 2`$.

## VI Universal Asymptotic Fluctuations

The fluctuations in the finite-data free energy, $`f_n=\mathrm{\Delta }F_n/k_BT`$, as measured by the variance $`\sigma _n`$ of $`_n`$ of (6), are of considerable interest because of their potential to provide parameter-free extrapolative estimates of $`f_{\mathrm{}}=\mathrm{\Delta }F/k_BT`$ ; see also . The variance is given by

$$\left(\frac{\sigma _n}{k_BT}\right)^2=\left(\frac{_n\mathrm{\Delta }F_n}{k_BT}\right)^2=_{cn^a}^{\mathrm{}}dy\rho _n(y)\left[\mathrm{log}(1+y/cn^a)\right]^2(f_nf)^2.$$
(68)
For $`n\mathrm{}`$, it was pointed out in  that the simple, linear relation between $`f_nf`$ and $`\sigma _n^2`$ was *independent* of the distribution of work values — that is, *universal*. Here, we sketch the derivation for the long-tailed case when the second Boltzmann moment diverges.

To calculate the asymptotic behavior of the fluctuations (68) note first that second term $`(f_nf)^2`$ is necessarily of higher order than $`f_nf`$. For the crucial integral of (68), one finds

$`{\displaystyle _{cn^a}^{\mathrm{}}}dy\rho _n(y)\left[\mathrm{log}(1+y/cn^a)\right]^2`$ $``$ $`{\displaystyle \frac{1}{n^{\alpha 1}}}{\displaystyle \frac{C_1^{\mathrm{}}}{c^\alpha }}I^{},`$ (69)
$`I_\alpha ^{}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}du{\displaystyle \frac{1}{u^{\alpha +1}}}\left[\mathrm{log}(1+u)\right]^2.`$ (70)
Comparing (68) - (70) with (67) and (66), we see that as $`n\mathrm{}`$

$$f_nf\frac{I_\alpha ^{}}{I_\alpha ^{}}\left(\frac{\sigma _n}{k_BT}\right)^2.$$
(71)
This is a linear relation that depends only on $`\alpha `$, via the ratio $`I_\alpha ^{}/I_\alpha ^{}`$, but is otherwise independent of the initital distribution of work values (or Boltzmann factors). In the limit $`\alpha 2`$, the ratio $`I_\alpha ^{}/I_\alpha ^{}`$ approaches $`1/2`$, which is the finite-Boltzmann-moment result reported in . Because numerical evaluation of the integral ratio is non-trivial we note that for $`\alpha =1.25`$, 1.5, 1.75, the corresponding values are $`I_\alpha ^{}/I_\alpha ^{}1.43`$, 0.81, 0.61.

Figure 3 illustrates the universal behavior for $`\alpha =1.5`$. Two integrable distributions were selected to ensure reliable computations. The “simple” or regulated-power-law distribution is defined by $`\rho _{rp}(z)=\alpha ^{}/(1+z)^{\alpha ^{}+1}`$, with $`\alpha ^{}=\alpha =1.5`$. The “power” distribution is given by $`\rho _p(z)=z_0/z^\alpha `$, with the choice $`z_0=10^4`$.

## VII Summary and Discussion

This report has expanded upon the brief discussion of Ref. , giving a general statistical theory describing the systematic error present in free-energy-difference $`\mathrm{\Delta }F`$ estimates based on a finite amount of data ($`N`$ work values, $`W`$). As in , our focus has been on the large-$`N`$ asymptotic behavior, motivated by the need to improve extrapolation procedures first explored in . However, beyond simply giving further details of the derivations of previous results, this report has made transparent the connection to general non-linear averages: the bounds of Sec. III, which generalize Jensen’s inequality, explicitly apply to a broad class of nonlinear computations in addition to $`\mathrm{\Delta }F`$ estimates; and, Sec. IV gives asymptotic expansions for geometric averages, such as the root-mean-square.

The universal, asymptotic relation (71) between the expected value of the biased $`\mathrm{\Delta }F`$ estimate based on $`N`$ work values ($`\mathrm{\Delta }F_N`$) and the fluctuation in these estimtates ($`\sigma _N`$) is one of the more striking results. We have shown here, in Sec. VI, that the relation is universal whether or not the second moment of the distribution of Boltzmann factors, $`\mathrm{exp}(W/k_BT)`$, is finite — that is, whether or not the central limit theorem applies. If not, the stable (Lévy) distributions come into play, and the relation between $`\mathrm{\Delta }F_N`$ and $`\sigma _N`$ depends only on the exponent of the limiting stable distribution.

We hope our results will have practical application in the extrapolation process outlined in , which suggested that dramatic increases in computational efficiency may be possible. In this context, examination of Padé approximants to the asymptotic series, which can be constructed to also exhibit suitable small-$`N`$ behavior, may prove fruitful. We believe, finally, that the statistical foundation laid in Ref.  and here provides a basis for the crucial but non-trivial task of simply *understanding* convergence in estimates of free energy differences and other non-linear averages.

###### Acknowledgements.

The authors have benefitted greatly from discussions and correspondence with Michael E. Fisher, Gerhard Hummer, Chris Jarzynski, David Kofke, Nandou Lu, Hagai Meirovitch, Hirsh Nanda, Lawrence Pratt, Mark Robbins, Attila Szabo, and David Zuckerman. Marty Ytreberg provided important comments in the final preparation of the manuscript. We gratefully acknowledge funding provided by the NIH (under grant GM54782 to T.B.W. and NRSA GM20394 to D.M.Z.), the Bard Foundation, and the Department of Physiology.

FIGURE CAPTIONS

Figure 1. Finite-sampling errors in $`\mathrm{\Delta }F`$ estimates based on (a) Gaussian-distributed work values and (b) work values generated in a molecular-mechanics solubility comparison between the fatty acids palmitate and stearate. The irregular, staircase-shaped plots are the running estimates based on $`N`$ work values, while the smooth curves depict the *average* running estimates $`\mathrm{\Delta }F_N`$ (5) which are independent of the order in which the work values were generated. The average running estimates are also rigorous upper bounds on $`\mathrm{\Delta }F`$. The standard deviation of the zero-mean Gaussian distribution in (a) is $`4k_BT`$, for which the true free energy difference is $`\mathrm{\Delta }F=\mathrm{\Delta }F_{\mathrm{}}=8k_BT`$. For the fatty acid solvation case, $`\mathrm{\Delta }F=\mathrm{\Delta }F_{\mathrm{}}13`$ kcal/mole; note that 1 kcal/mole = 1.7 $`k_BT`$.

Figure 2. Finite-sampling error for Gaussian-distributed work values. The expected value of the dimensionless finite-sampling inaccuracy, $`(\mathrm{\Delta }F_n\mathrm{\Delta }F)/k_BT`$ for $`n`$ data points is plotted as a function of $`1/n`$. From top to bottom, the data sets represent numerical values of the error for Gaussian distributions of work values with standard deviations, $`\sigma _w/k_BT`$ of 3, 2, 1.5, and 1. The lines (dashed for $`\sigma _w/k_BT=1.5`$, solid for $`\sigma _w/k_BT=1`$) depict the asymptotic linear behavior for the two smallest widths.

Figure 3. The universal $`n\mathrm{}`$ relation between $`\mathrm{\Delta }F_n\mathrm{\Delta }F`$ and its fluctuation $`\sigma _n`$ for the long-tailed case when the Lévy index is $`\alpha =1.5`$. The solid line depicts the universal slope $`I_\alpha ^{}/I_\alpha ^{}0.815`$ for $`\alpha =1.5`$, as given in (71) and the succeeding text. The data for the “power” and “simple” distributions, described in the text, are each shown for $`\alpha ^{}=\alpha =1.5`$.