% PRL look and style (easy on the eyes)
%\documentclass[aps,prl,twocolumn,superscriptaddress]{revtex4}
% Two-column style (for submission/review/editing)
\documentclass[aps,preprint,prl,superscriptaddress,showpacs]{revtex4}

\usepackage{palatino}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{dcolumn}
\usepackage{boxedminipage}
\usepackage{verbatim}
\usepackage[colorlinks=true,citecolor=blue,linkcolor=blue]{hyperref}

% The figures are in a figures/ subdirectory.
\graphicspath{{figures/}}

%\bibliographystyle{apsrevlong}
%\bibliographystyle{apsrev}

% italicized boldface for math (e.g. vectors)
\newcommand{\bfv}[1]{{\mbox{\boldmath{$#1$}}}}
% non-italicized boldface for math (e.g. matrices)
\newcommand{\bfm}[1]{{\bf #1}}          

%\newcommand{\bfm}[1]{{\mbox{\boldmath{$#1$}}}}
%\newcommand{\bfm}[1]{{\bf #1}}
\newcommand{\expect}[1]{\left \langle #1 \right \rangle}                % <.> for denoting expectations over realizations of an experiment or thermal averages

% vectors
\newcommand{\x}{\bfv{x}}
\newcommand{\y}{\bfv{y}}
\newcommand{\f}{\bfv{f}}

\newcommand{\bfc}{\bfm{c}}
\newcommand{\hatf}{\hat{f}}

\newcommand{\bTheta}{\bfm{\Theta}}
\newcommand{\btheta}{\bfm{\theta}}
\newcommand{\bhatf}{\bfm{\hat{f}}}
\newcommand{\Cov}[1] {\mathrm{cov}\left( #1 \right)}
\newcommand{\Var}[1] {\mathrm{var}\left( #1 \right)}
\newcommand{\Ept}[1] {{\mathrm E}\left[ #1 \right]}
\newcommand{\dhdl}{\frac{dH}{d\lambda}}
\newcommand{\Eptk}[2] {{\mathrm E}_{#1}\left[ #2\right]}
\newcommand{\T}{\mathrm{T}}                                % T used in matrix transpose

\newcommand{\pymbar}{{\em pymbar}}
%% DOCUMENT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

%% TITLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Applications and advantages of the multistate Bennet acceptance ratio estimato}

\author{John D. Chodera}
 \email{jchodera@berleley.edu}
 \affiliation{Center for Quantitative Biololgy, University of California, Berkeley, CA}
\author{Michael R. Shirts}
  \email{michael.shirts@virgina.edu}
 \affiliation{Department of Chemical Engineering, University of Virginia, Charlottesville, VA 22094-0471}

\date{\today}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ABSTRACT/pacs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}

Stuff here.

\end{abstract}

\pacs{02.30.Cj,05.70.Ce,82.20.Wt}
\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INTRODUCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

A recurring challenge in statistical physics, computational chemistry,
and single molecule experiments is the collection of a sufficient
amount of data to estimate physical quantities of interest to adequate
precision.  In computer simulations of physical or chemical models,
such quantities include potentials of mean force, phase coexistence
curves, fluctuation or temperature-dependent properties, and free
energy differences.  In single-molecule experiments, these quantities
might include potentials of mean force along a pulling coordinate or
the distance between fluorescence probes during resonant energy
transfer.  For these problems, collection of sufficient statistics for
a reliable estimate often requires multiple simulations at different
thermodynamic states~\cite{MLEmeasure_foot1} or measurements performed
under different applied fields.  In computer simulations, multi-state
techniques such as umbrella
sampling~\cite{torrie-valleau:j-comput-phys:1977:umbrella-sampling},
simulated~\cite{marinari-parisi:europhys-lett:1992:simulated-tempering}
and parallel
tempering~\cite{hukushimi-nemoto:j-phys-soc-jpn:1996:parallel-tempering}
and the use of alchemical intermediates in free energy calculations
can greatly aid convergence.

Even with these methods, it may require a large quantity of data to produce estimates with the desired precision. 
Computing the most precise estimate possible from the available data can therefore be critical in allowing these quantities to be estimated with reasonable computational or experimental effort. 
While the choice of thermodynamic states to sample can also greatly affect efficiency, we focus here on only the problem of statistically efficient estimation given samples from predetermined states.
Early methods for computing free energy differences~\cite{zwanzig:1954a,widom:1963} or equilibrium expectations~\cite{torrie-valleau:j-comput-phys:1977:umbrella-sampling} relied upon one-sided exponential averaging (EXP), which is formally exact but does not make the most efficient use of data when samples from more than one state are available~\cite{shirts:2005a}.
Subsequently, the Bennett acceptance ratio method (BAR)~\cite{bennett:1976,shirts:2003a} greatly improved upon EXP for the computation of free energy differences, producing statistically optimal estimates of free energy differences when two states are sampled~\cite{shirts:2003a} and yielding estimates that can be more than an order of magnitude more precise~\cite{shirts:2005a}.
More recently, multiple histogram reweighting methods~\cite{ferrenberg:1989a,kumar:1992a} were proposed as a way to incorporate data from multiple states to produce superior estimates of free energy differences and equilibrium expectations for arbitrary thermodynamic states, including states not sampled.

While multiple histogram techniques --- most notably, the weighted histogram analysis method (WHAM)~\cite{kumar:1992a} --- can produce statistically optimal estimates of the discretized densities of states~\cite{ferrenberg:1989a} or histogram occupation probabilities~\cite{bartels:1997a}, they have several limitations for the treatment of continuous systems.
First, the reliance on histograms of width sufficient to contain many samples introduces a bias that can be substantial and often difficult to assess~\cite{kobrak:j-comput-chem:2003:histogram-error}.  
Second, unlike BAR, there are no direct expressions to estimate the statistical uncertainty in free energy differences or expectations obtained from WHAM.
Third, application of WHAM to a samples collected with a biasing potential that is not trivially scaled by a linear field parameter (such as inverse temperature) requires a number of bins that grows exponentially in the number of states, making it computationally intractable for even modest numbers of states.
While more recent maximum likelihood~\cite{bartels:1997a} and Bayesian formulations~\cite{gallicchio:2005a} mitigate the memory requirements, they do not remove the histogram bias effects, and introduce a costly Markov chain Monte Carlo sampling procedure to estimate uncertainties~\cite{gallicchio:2005a,park:pre:2006:abwham}.

Here, we use recent results from the field of statistical inference~\cite{vardi:1985,gill:1988,kong:2003,tan:2004} to construct a statistically optimal estimator for computing free energy differences and equilibrium expectations at arbitrary thermodynamic states, using equilibrium samples from multiple thermodynamic states.
The resulting estimator, termed the multistate Bennett acceptance ratio (MBAR) estimator as it reduces to Bennett's method when only two states are considered~\cite{bennett:1976}, is equivalent to WHAM in the limit that histogram bin widths are shrunk to zero but is derived without the need to invoke histograms.  
Unlike WHAM, this estimator provides a direct assessment of uncertainties, critical in making comparison between experiment and theory, and the computational expense of computing the estimator remains modest across a wider variety of applications.
Furthermore, it can easily be applied to data sampled from non-Boltzmann sampling schemes, or to the analysis of single-molecule experiments in cases where an external bias potential is applied.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXTENDED BRIDGE SAMPLING ESTIMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Extended bridge sampling estimators}

Suppose we obtain $N_i$ uncorrelated equilibrium samples~\cite{MLEmeasure_foot2} from each of $K$ thermodynamics states within the same ensemble (such as NVT, NPT, or $\bfm{\mu}$VT).
Each state is characterized by a specified combination of inverse temperature, potential energy function, pressure, and/or chemical potential(s), depending upon the ensemble.
We define the \emph{reduced potential function} $u_i(\x)$ for state $i$ to be
\begin{eqnarray}
u_i(\x) &=& \beta_i [ U_i(\x) + p_i V(\x) + \bfv{\mu}_i^\T \bfv{n}(\x) ] \label{equation:reduced-potential}
\end{eqnarray}
where $\x \in \bfm{\Gamma}$ denotes the configuration of the system within a configuration space $\bfm{\Gamma}$, with volume $V(\x)$ (in the case of a constant pressure ensemble) and $\bfv{n}(\x)$ the number of molecules of each of $M$ components of the system (in the case of a (semi)grand ensemble).  
For each state $i$, $\beta_i$ denotes the inverse temperature, $U_i(\x)$ the potential energy function (which may include biasing weights), $p_i$ the external pressure, and $\bfv{\mu}_i$ the vector of chemical potentials of the $M$ system components.  

Configurations $\{\x_{in}\}_{n=1}^{N_i}$ from state $i$ are sampled from the probability distribution 
\begin{eqnarray}
p_i(\x) = c_i^{-1} q_i(\x) \: &;& \: c_i = \int_{\bfm{\Gamma}} d\x \, q_i(\x)
\end{eqnarray}
where $q_i(\x)$ is here nonnegative and represents an unnormalized density function, and $c_i$ is the (generally unknown) normalization constant (known in statistical mechanics as the \emph{partition function}).
In samples obtained from standard Metropolis Monte Carlo or molecular dynamics simulations or from experiment, this unnormalized density is simply the Boltzmann weight $q_i(\x) = \exp[-u_i(\x)]$ but may in general differ in simulations employing non-Boltzmann weights, such as multicanonical simulations~\cite{mezei:j-comp-phys:1987:muca} and those using Tsallis statistics~\cite{tsallis:j-stat-phys:1988:tsallis-statistics}.

We wish to produce an estimator for the difference in dimensionless free energies
\begin{eqnarray}
\Delta f_{ij} &\equiv& f_j - f_i = - \ln \frac{c_j}{c_i} = - \ln \frac{\int_\bfm{\Gamma} d\x \, q_j(\x)}{\int_\bfm{\Gamma} d\x \, q_i(\x)} \label{equation:dimensionless-free-energy}
\end{eqnarray}
where the $f_i$ are related to the dimensional free energies $F_i$ by $f_i = \beta_i F_i$, and also the equilibrium expectations
\begin{eqnarray}
\expect{A}_i &\equiv& \int_{\bfm{\Gamma}} d\x \, p_i(\x) \, A(\x) = \frac{\int_\bfm{\Gamma} d\x \, A(\x) \, q_i(\x)}{\int_\bfm{\Gamma} d\x \, q_i(\x)} .\label{equation:equilibrium-expectation}
\end{eqnarray}
These expectations can be computed as ratios of the normalization constants if we define new functions $q(\x) = A(\x) q_i(\x)$, where the $q(\x)$ no longer need be nonnegative for states from which no samples are collected~\cite{doss:2003}.

To construct an estimator for these ratios of normalization constants, we first note the identity
\begin{eqnarray}
c_i \expect{\alpha_{ij} q_j}_i = \int_{\bfm{\Gamma}} d\x \, q_i(\x) \, \alpha_{ij}(\x) \, q_j(\x)  =  c_j \expect{\alpha_{ij} q_i}_j
\label{equation:Bennett_starting}
\end{eqnarray}
which holds for arbitrary choice of functions $\alpha_{ij}(\x)$, provided the $c_i$ are nonzero.

Using this relation, summing over the index $j$, and substituting the empirical estimator $N_i^{-1} \sum_{n=1}^{N_i} g(\x_{in})$ for the expectations $\expect{g}_i$, we obtain a set of $K$ estimating equations
\begin{eqnarray}
\sum_{j=1}^K \frac{\hat{c}_i}{N_i} \sum_{n=1}^{N_i} \alpha_{ij} q_j(\x_{in}) &=& \sum_{j=1}^K \frac{\hat{c}_j}{N_j} \sum_{n=1}^{N_j} \alpha_{ij} q_i(\x_{jn}) \label{equation:extended-bridge-sampling-estimator}
\end{eqnarray}
for $i = 1,2,\ldots,K$, where solution of the set of equations for the $\hat{c}_i$ yields an estimate of the $c_i$ from the sampled data determined up to a scalar multiplier.

Eq.~\ref{equation:extended-bridge-sampling-estimator} defines a family of asymptotically unbiased estimators parameterized by the choice of functions $\alpha_{ij}(\x)$, known in the statistics literature as \emph{extended bridge sampling} estimators~\cite{tan:2004}.
By making the choice
\begin{eqnarray}
\alpha_{ij}(\x) &=& N_j \, \hat{c}_j^{-1} / \sum\limits_{k=1}^K N_k \, \hat{c}_k^{-1} \, q_k(\x) \:\:, \label{equation:optimal-alpha}
\end{eqnarray}
we obtain an estimator than has been proven to be optimal, in the sense that it has the lowest variance of for all choices of $\alpha_{ij}(\x)$~\cite{tan:2004}.  
This estimator is also asymptotically unbiased and guaranteed to have a unique solution (up to a multiplicative constant)~\cite{tan:2004}, and can also been derived from maximum-likelihood measure theoretic methods~\cite{kong:2003} or cast as a reverse logistic regression problem~\cite{geyer:1994,kong:2003}.

While a closed-form expression for the $\hat{c}_i$ cannot be obtained from Eqs.~\ref{equation:extended-bridge-sampling-estimator} and \ref{equation:optimal-alpha}, the $\hat{c}_i$ can nevertheless be easily computed by any suitable method for solving systems of coupled nonlinear equations.
A simple self-consistent iteration method and an efficient Newton-Raphson solver are described in the Supplementary Information, and Python code implementing these methods can be found at https://simtk.org/home/pymbar.

In the large sample limit, the error in the ratios $\hat{c}_i/\hat{c}_j$ will be normally distributed~\cite{tan:2004}, and the asymptotic covariance matrix $\bfm{\Theta}_{ij} = \Cov{\theta_i, \theta_j}$, where $\theta_i \equiv \ln c_i$, can be estimated by~\cite{kong:2003}
\begin{eqnarray}
\hat{\bfm{\Theta}} &=& \bfm{W}^\T (\bfm{I}_N - \bfm{W} \bfm{N} \bfm{W}^\T)^{+} \bfm{W} \label{equation:asymptotic-covariance}
\end{eqnarray}
where $\bfm{I}_N$ is the $N \times N$ identity matrix  (with $N = \sum_{i=1}^K N_i$ the total number of samples), and $\bfm{N} = \mathrm{diag}(N_1, N_2,\ldots,N_K)$.
The superscript $^+$ denotes a suitable generalized inverse, such as the Moore-Penrose pseudoinverse, since the quantity in parentheses will be rank-deficient.
$\bfm{W}$ denotes the $N \times K$ matrix of weights
\begin{eqnarray}
W_{ni} = \hat{c}_i^{-1} \, q_i(\x_n) / \sum\limits_{k=1}^K N_k \, \hat{c}_k^{-1} \,\ q_k(\x_n) \label{equation:p_weights} .
\end{eqnarray}
The samples are now indexed by a single index $n = 1,\ldots,N$, as the association of which samples $\x_n$ came from which distribution $p_i(\x)$ is no longer relevant. 
The computational cost of computing the pseudoinverse of an $N \times N$ matrix in computing $\hat{\bfm{\Theta}}$ can be reduced to that of computing the eigenvalue decomposition of a $K \times K$ matrix (see Supplementary Information).

The covariance of estimates of arbitrary functions $\phi(\theta_1,\ldots,\theta_K)$ and $\psi(\theta_1,\ldots,\theta_K)$ of the normalization constants $c_i$ can be estimated from $\hat{\bfm{\Theta}}$ by the expansion
\begin{eqnarray}
\Cov{\hat{\phi},\hat{\psi}} &=& \sum_{i,j = 1}^K \frac{\partial \phi}{\partial \theta_i} \,\hat{\Theta}_{ij} \frac{\partial \psi}{\partial \theta_j} \,.\label{equation:asymptotic-covariance-transformation}
\end{eqnarray}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ESTIMATOR FOR FREE ENERGY DIFFERENCES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Free energies}

When configurations are sampled with Boltzmann statistics, Eqs.~\ref{equation:extended-bridge-sampling-estimator} and \ref{equation:optimal-alpha} produce the following estimating equations for the dimensionless free energies
\begin{eqnarray}
\hat f_i &=& - \ln \sum_{j=1}^K \sum_{n=1}^{N_j} \frac{\exp[-u_i(\x_{jn})]}{\sum\limits_{k=1}^K N_k \, \exp[\hat{f}_k - u_k(\x_{jn})]} \label{equation:estimator-of-free-energies}
\end{eqnarray}
which must be solved self-consistency for the $\hat{f}_i$.
Again, because the normalization constants are only determined up to a multiplicative constant, the estimated free energies $\hat{f}_i$ are determined uniquely only up to an additive constant, so only differences $\Delta \hat{f}_{ij} = \hat{f}_j - \hat{f}_i$ will be meaningful.

The uncertainty in the estimated free energy difference can be computed from Eqs.~\ref{equation:asymptotic-covariance} and \ref{equation:asymptotic-covariance-transformation} as
\begin{eqnarray}
\delta^2 \Delta \hat{f}_{ij} &\equiv& \Cov{-\ln \hat{c}_j/\hat{c}_i, -\ln \hat{c}_j/\hat{c}_i} = \hat{\Theta}_{ii} - 2 \hat{\Theta}_{ij} + \hat{\Theta}_{jj} \nonumber
\end{eqnarray}
Free energy differences and uncertainties between states not sampled are easily estimated by augmenting the set of states with additional reduced potentials $u_i(\x)$ with the number of samples $N_i = 0$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ESTIMATOR FOR EXPECTATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Expectations}

The equilibrium expectation of some mechanical observable $A(\x)$ that depends only on configuration $\x$ (and not momentum) is given by Eq.~\ref{equation:equilibrium-expectation}, and can be computed as a ratio of normalization constants $c_A/c_a$ by defining two additional functions
\begin{eqnarray}
q_A(\x) = A(\x) \, q(\x) \:&;&\: q_a(\x) = q(\x) \nonumber 
\end{eqnarray}
where $q(\x) = \exp[-u(\x)]$ if the expectation with respect to the Boltzmann weight is desired~\cite{doss:2003}.
While $q_A(\x)$ may no longer be strictly nonnegative, we may still make use of the extended bridge sampling estimator (Eq.~\ref{equation:extended-bridge-sampling-estimator}) since $N_A = N_a = 0$.

If we augment the matrix $\bfm{W}$ (Eq.~\ref{equation:p_weights}) with the columns $W_{nA}$ and $W_{na}$ corresponding to $q_A(\x)$ and $q_a(\x)$, respectively, we can write the estimator as
\begin{eqnarray}
\hat{A} = \hat{c}_A / \hat{c}_a = \sum\limits_{n=1}^N W_{nA} / \sum\limits_{n=1}^N W_{na} = \sum_{n=1}^N W_{na} \, A(\x_{n}) . \label{equation:estimator-for-expectations}
\end{eqnarray}
We note that, if the dimensionless free energies $\hat{f}_i$ have already been determined, computation of $\hat{A}$ for any $A(\x)$ and any $q(\x)$ does not require additional solution of self-consistent equations.

The uncertainty in $\hat{A}$ is given by
\begin{eqnarray}
\delta^2 \hat{A} \equiv \Cov{\hat{c}_A/\hat{c}_a,\hat{c}_A/\hat{c}_a} = \hat{A}^2 (\hat{\Theta}_{AA} + \hat{\Theta}_{aa} -  2 \, \hat{\Theta}_{Aa}) \nonumber
\end{eqnarray}
where the covariance matrix $\hat{\bfm{\Theta}}$ is computed from the augmented $\bfm{W}$.
Covariances between estimates of $\expect{A}$ at different thermodynamic states, or between two observables $\expect{A}$ and $\expect{B}$, can also be constructed by adding the appropriate columns to the covariance matrix and applying Eq.~\ref{equation:asymptotic-covariance-transformation} to estimate the desired uncertainty.  
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISCUSSION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Discussion}

The MBAR estimator presented here provides a rapid and robust way to extract estimates of free energy differences and equilibrium expectations from multiple equilibrium samples of different thermodynamic states in a statistically optimal way.
As the estimator is asymptotically efficient among a wide class of ``bridge sampling'' estimators~\cite{tan:2004}, which includes EXP and BAR as members, the resulting estimates from MBAR will have the lowest (or equal) variance in the large sample limit.

While multiple histogram techniques~\cite{kumar:1992a,ferrenberg:1989a,bartels:1997a,gallicchio:2005a} have been widely used for combining data from multistate simulations, the MBAR estimator supplants, or at least complements, these methods.
Perhaps most importantly, it provides a reliable and inexpensive method for estimating the uncertainties in the resulting estimates and their correlations, which are critical for propagating uncertainties to quantities of interest.  
Additionally, the elimination of histograms avoids both the bias arising from discretization of continuous energies as well as the computational overhead of constructing and storing high dimensional histograms.  

In this framework, multiple histogram reweighting methods such as WHAM can be understood as a histogram kernel density estimator approximation to MBAR.
In some applications, histograms can reduce the computational expense required for solving the estimating equations (Eq.~\ref{equation:extended-bridge-sampling-estimator}) at the expense of introducing bias. 
When samples are distributed according to the Boltzmann weight, the estimator for the free energies (Eq.~\ref{equation:estimator-of-free-energies}) is precisely Eq.~21 of~\cite{kumar:1992a} or Eq.~15 of~\cite{souaille:2001a}, in both cases presented as a reduction of the histogram bin width to zero in the standard WHAM equations (Eqs.~19--20 of \cite{kumar:1992a}).
While the validity of this limit is dubious --- the derivations in these references rely upon an estimate of the uncertainty in each histogram count which cannot be correct when the bins are nearly empty --- the derivation of this equation from the extended bridge sampling estimator demonstrates for the first time that these equations are, in fact, asymptotically unbiased estimators of the true free energy differences.  

The MBAR estimator also can be considered a multistate generalization of Bennett's acceptance ratio estimator (BAR)~\cite{bennett:1976}.  
In deriving BAR, Bennett constructed an estimator from Eq.~\ref{equation:Bennett_starting} directly, determining the single $\alpha(\x)$ which minimized the variance of the estimator of the free energy difference between only two states.
In deriving MBAR, summing over all states $j$ and determining the functions $\alpha_{ij}(\x)$ that minimizes the covariance matrix of the estimator for ratios of normalization constants produces an optimal estimator for the multistate case.  
A proof of the equivalence of MBAR and BAR for two states is in the Supplementary Information.

BAR and a recent pairwise multistate generalization (which we shall refer to as PBAR)~\cite{maragakis:prl:2006} differ from MBAR in that they can also be applied to \emph{nonequilibrium} work measurements between pairs of states, in addition to equilibrium reduced potential differences (\emph{instantaneous} work measurements).  
However, PBAR constructs a total likelihood function from products of likelihood functions connecting pairs of states, assuming independence of all work measurements.  
For equilibrium samples, this means that a sampled configuration $\x_{n}$ from a state $i$ can only be used to provide information about the instantaneous work required to switch to a \emph{single} other state $j$ for use in the PBAR estimator, whereas in MBAR, each sampled $\x_{n}$ can be used to provide information about \emph{all} states.  
As a result, MBAR should require significant fewer samples from each state to produce an estimate of equivalent precision.


Finally, we note the MBAR estimator is not limited in application to data produced from simulation --- it can also be applied to combine data from multiple equilibrium experiments in the presence of externally-applied fields.
To illustrate, we used MBAR to estiamate potential of mean force (PMF) of a DNA hairpin attached by dsDNA linkers to glass beads along the distance between the beads.
A set of 16 equilibrium trajectories of 50 000 samples each under a variety of constant force loads (corresponding to a linear external potential along the extension coordinate) for the hairpin system 20R55/4T collected by an optical double trap experiment~\cite{woodside:prl:2005:optical-force-clamp,woodside:pnas:2006:dna-hairpin-equilibrium} were analyzed to produce an optimal estimate of the PMF under a force load of 14.19 pN (Supplementary Information, Fig.~1).
The analysis took 18 s on a standard 2.16 GHz Intel Core 2 Duo MacBook Pro, and the resulting error bars are more than an order of magnitude smaller than those derived from the single trajectory at this force load in the poorly sampled region of the PMF.

\input alchemical.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ACKNOWLEDGMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Acknowledgments}

We are indebted to Michael Woodside for providing us with detailed datasets from Ref.~\cite{woodside:pnas:2006:dna-hairpin-equilibrium} and helpful comments.
We thank Evangelos A.~Coutsias, Gavin E.~Crooks, Fernando A.~Escobedo, Edward H.~Feng, Andrew Gelman, Andrew I.~Jewett, Libusha Kelly, Jun S.~Liu, David D.~L.~Minh,  David L.~Mobley, Frank M.~No\'{e}, Vijay S.~Pande, Sanghyun Park, M.~Scott Shell, Zhiqiang Tan, Matthew A.~Wyczalkowski, and the anonymous reviewers for enlightening discussions and constructive comments on this manuscript. 
JDC gratefully acknowledges support from Ken A.~Dill through NIH grant GM34993 and Vijay S.~Pande through an NSF grant for Cyberinfrastructure (NSF CHE-0535616), and MRS support from Richard A.~Friesner and an NIH NRSA Fellowship.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BIBLIOGRAPHY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{prsty} 

\begin{thebibliography}{10}

\bibitem{MLEmeasure_foot1}
Here, a \emph{thermodynamic state} is defined by a combination of potential
  energy function (including any biasing potentials) and external thermodynamic
  parameters, such as temperature, pressure, and chemical potential, all within
  the same thermodynamic ensemble (e.g. canonical, isothermal-isobaric,
  (semi)grand canonical).

\bibitem{torrie-valleau:j-comput-phys:1977:umbrella-sampling}
G.~M. Torrie and J.~P. Valleau, J. Comput. Phys. {\bf 23},  187  (1977).

\bibitem{marinari-parisi:europhys-lett:1992:simulated-tempering}
E. Marinari and G. Parisi, Europhys. Lett. {\bf 19},  451  (1992).

\bibitem{hukushimi-nemoto:j-phys-soc-jpn:1996:parallel-tempering}
K. Hukushimi and K. Nemoto, J. Phys. Soc. Jpn. {\bf 65},  1604  (1996).

\bibitem{zwanzig:1954a}
R.~W. Zwanzig, J. Chem. Phys. {\bf 22},  1420  (1954).

\bibitem{widom:1963}
B. Widom, J. Chem. Phys {\bf 39},  2808  (1963).

\bibitem{shirts:2005a}
M.~R. Shirts and V.~S. Pande, J. Chem. Phys. {\bf 122},  144107  (2005).

\bibitem{bennett:1976}
C.~H. Bennett, J. Comp. Phys. {\bf 22},  245  (1976).

\bibitem{shirts:2003a}
M.~R. Shirts, E. Bair, G. Hooker, and V.~S. Pande, Phys. Rev. Lett. {\bf 91},
  140601  (2003).

\bibitem{ferrenberg:1989a}
A.~M. Ferrenberg and R.~H. Swendsen, Phys. Rev. Lett. {\bf 63},  1195  (1989).

\bibitem{kumar:1992a}
S. Kumar, D. Bouzida, R.~H. Swendsen, P.~A. Kollman, and J.~M. Rosenberg, J.
  Comput. Chem. {\bf 13},  1011  (1992).

\bibitem{bartels:1997a}
C. Bartels and M. Karplus, J. Comput. Chem. {\bf 18},  1450  (1997).

\bibitem{kobrak:j-comput-chem:2003:histogram-error}
M.~N. Kobrak, J. Comput. Chem. {\bf 24},  1437  (2003).

\bibitem{gallicchio:2005a}
E. Gallicchio, M. Andrec, A.~K. Felts, and R.~M. Levy, J. Phys. Chem. B {\bf
  109},  6722  (2005).

\bibitem{park:pre:2006:abwham}
S. Park, D.~L. Ensign, and V.~S. Pande, Phys. Rev. E {\bf 74},  066703  (2006).

\bibitem{vardi:1985}
Y. Vardi, Ann. Stat. {\bf 13},  178  (1985).

\bibitem{gill:1988}
R.~D. Gill, Y. Vardi, and J.~A. Wellner, Ann. Stat. {\bf 16},  1069  (1988).

\bibitem{kong:2003}
A. Kong, P. McCullagh, X.-L. Meng, D. Nicolae, and Z. Tan, J. Royal Stat. Soc.
  B. {\bf 65},  585  (2003).

\bibitem{tan:2004}
Z. Tan, J. Am. Stat. Assoc. {\bf 99},  1027  (2004).

\bibitem{MLEmeasure_foot2}
A set of uncorrelated configurations can be obtained from a correlated time
  series, such as is generated by a molecular dynamics or Metropolis Monte
  Carlo simulation, by subsampling the timeseries (see Supplementary Information).

\bibitem{mezei:j-comp-phys:1987:muca}
E. Mezei., J. Comp. Phys. {\bf 68},  237  (1987).

\bibitem{tsallis:j-stat-phys:1988:tsallis-statistics}
C. Tsallis, J. Stat. Phys. {\bf 52},  479  (1988).

\bibitem{doss:2003}
Honi Doss makes this suggestion in the conference discussion
  of~\protect\cite{kong:2003}.

\bibitem{geyer:1994}
C.~J. Geyer, Technical Report No.~568, School of Statistics, University of
  Minnesota, Minneapolis, Minnesota (unpublished).

\bibitem{souaille:2001a}
M. Souaille and B. Roux, Comp. Phys. Commun. {\bf 135},  40  (2001).

\bibitem{maragakis:prl:2006}
P. Maragakis, M. Spichty, and M. Karplus, Phys. Rev. Lett. {\bf 96},  100602
  (2006).

\bibitem{woodside:prl:2005:optical-force-clamp}
W.~J. Greenleaf, M.~T. Woodside, E.~A. Abbondanzieri, and S.~M. Block, Phys.
  Rev. Lett. {\bf 95},  208102  (2005).

\bibitem{woodside:pnas:2006:dna-hairpin-equilibrium}
M.~T. Woodside, W.~M. Behnke-Parks, K. Larizadeh, K. Travers, D. Herschlag, and
  S.~M. Block, Proc. Nat. Acad. Sci. USA {\bf 103},  6190  (2006).

\end{thebibliography}

\end{document}
