\section{Validation of MBAR formalism}

The main theoretical claims of MBAR are that 1) in the limit of large
numbers of samples, the errors in free energies or ensemble averages
are normally distributed with zero bias and with variance given by the
the MBAR analytical estimate of variance and 2) MBAR is the
reweighting method with the lowest asymptotic variance.

To demonstrate these properties of MBAR and its asymptotic variance
and to validate the implementation of MBAR in {\em pymbar}, we analyze
the case of samples taken from systems of offset harmonic oscillators.
Systems of offset harmonic oscillators are simple enough to analyze
analytically, and flexible enough to test a range of parameters and
situations.  Of course, they are simplified system, and will not
entirely reflect the behavior of MBAR or other free energy and
expectation estimators in other situations.

We first examine the claim that the errors of free energies and
ensembles averages are normally distributed with zero bias in the
limit of large numbers of samples, and that the variance equal to the
analytical estimate of the variance.  We can test the validity of this
by repeating the simulation, drawing new samples from each
thermodynamic state $N$ times.  For each repetition of the experiment,
we can calculate the error $\delta_i$ in the calculated free energy
difference as the difference between the MBAR estimate of the free
energy or expectation, and the analytical estimate of the free
energy. We then divide by the MBAR estimate of the variance from this
repetition of the experiment $\hat{\sigma}_i$ to obtain the estimated
number of standard deviations of the estimate $s_i$ =
$\delta_i/\hat{\sigma}_i$.  If the distribution of error is normal,
then the cumulative distribution function (CDF) of $s$ should simply
be the normal distribution's CDF, the error function.  Specifically,
$s$ should be less than 1 $\hat{\sigma}$ 68\% of the time, less than 2
$\hat{\sigma}$ 95\% if the time, and so forth.

If the distribution of errors is Gaussian, but variance estimate is an
{\em underestimate} of the true variance, then the CDF of $s(X)$ will
be uniformly lower than the normal distribution curve.  If the
distribution of errors is Gaussian, then the variance estimate is an
{\em overestimate} of the true variance, then the CDF of $s(X)$ will
be uniformly higher than normal distribution curve.  If the
distribution is non-Gaussian, then the deviations will be more
complicated, and could involve underestimates in some parts of the
curve, and overestimates in other parts.  In general, we are most
concerned with long tail behavior, as it means that errors have some
nonneligible chance to be much greater than the variance estimate.
This shows up as the CDF of $s$ undershooting of the normal curve in
the large $\sigma$ regime.

For each of the simulations below, we set $\beta=1$, and take as the
reference the state at $x=0$. For reference he reduced free energy of
the oscillators will be $f_i = -\ln\sqrt{2\pi/\beta k_i}$, the average
reduced potential energy will be $\expect{u_i} = (2\beta)^{-1}$, the
average position will be $x_i$, and the average square position will
be $(\beta k_i)^{-1}+ x_i^2$.

\subsubsection{Equal spring constants}

We first examine a system of seven harmonic oscillators, with
equilibrium positions $x^0_i = i-1$, and spring constants $k_i = 16$,
with $\beta = 1$, as a function of $N$, the number of samples from
each state, equal for all oscillators.  We choose $k=16$, since
$\sigma = (\beta k)^{-1}$ so that each distribution is $2\sigma$ from
the midpoint between neighboring oscillators.  This arangements
represents small but nonnegligible phase space overlap between
neighbors; on average, only 5\% of the displacements will be closer to
the equilibrium value of their neighbor than they will be to their own
equilibrium value, and $\int_{-\infty}^{\infty} p_i(x) p(x) dx \approx
0.10$.  We examine both the free energies and the
$\expect{x^2}$. 

In Fig.~\ref{fig:cdfequalspace_N50}, we graph the CDF's of
$s(f_i-f_1)$ for $i=2,\ldots,7$ for $N=50$. We perform 10,000 replicates
of the experiment, so the error bars are essentially negligible, and
graph the uncertainties of

In Fig.~\ref{fig:cdfequalspace_bynumber}, we graph the value
of $s(f_i-f_1)$ at 2 $\sigma$ as a function of $N$.

\subsubsection{Unequal spring constants}

This performace is encouraging, but there are a number of symmetries
that might result in cancelling errors when we have identical harmonic
oscillators.  We now look at the peformance with unequal spring
constants.  We use the same spacing as before, $x^0_i = i-1$, but now
choose $k_i = 4,8,16,32,16,8,4$.

\subsubsection{Unsampled states}

One of the main features of MBAR is the ability to compute free
energies and ensemble averages for states that are not sampled.  To
test the performance of the estimators in this case, We use the same
unequal spring constants as before, but fix the total number sampled
at each state to $N=50$, and examine both the case that $N_4$ (the
largest spring constant) is unsampled, and that $N_7$ (the last state)
is unsampled, and look at the free energy estimates and expectation
values.

\subsubsection{Other observables}

An important part of the MBAR method is allowing the computation of
arbitrary observables using simulations at any thermodynamic state
using the simulations collected at the observable states in a
statistically efficient way.  Additionally, the uncertanties in these
observables can be computed, inplicitly assuming the uncertainties in
the free energies.

We look at the distribution of 4 observables for our system of 7
oscillators, each of which can be computed analytically: the average
positions of each oscillator, $\expect{x}_i = x^0_i$, the means square
deviation of each oscillator $\expect{(x-x^0_i)^2} = \sigma_i$, the
average potential energy $\expect{v}_i = (2\beta)^{-1}$, and the
average position squared, $\expect{x^2}_i = (\beta k_i)^{-1} +
(x^0_i)^2$.

As an example, we show the 
\begin{table}
\end{table}

In figure~\cite{???}


\subsubsection{Entropy-enthalpy decomposition}

We can also use the normal oscillators as a test of other
thermodynamically important estimators. A challenging computation of
importance in many types of molecular simulation is the decomposition
of free energy into entropic and energetic terms.  The reduced free
energy $f_i$, average reduced energy $u_i = \expect{v_i}_i$, and
reduced entropy $s_i = S/k_B$ obey the Helmholtz relationship:
\begin{eqnarray}
f_i &=& u_i - s_i \nonumber
\end{eqnarray}
A straightforward estimator for the reduced entropy $s_i$ is therefore given by
\begin{eqnarray}
\hat{s}_i &=& - \hat{f}_i + \hat{u}_i = -\hat{f}_i + \sum_{n=1}^N w_{ni} v_i(x_n) \label{equation:simple-entropy-estimator}
\end{eqnarray}

We take the case where $k_i = 0.2$.

For each experiment, 50 
independent and identically distributed
(\emph{i.i.d.}) 

samples were collected from each of $K = 7$
thermodynamic states, differing only in x0, where x0 was one of $\{0,
0.3, 0.6, 0.9, 1.2, 1.5, 1.8\}$.  1000 independent trials of this
experiment were run, and statistics were collected on the error in the
estimated entropy difference between state 1 ($x_0 = 0$) and state 7
($x_0 = 1.8$), where the true thermodynamic differences are
identically zero ($\Delta s_{17} = \Delta u_{17} = \Delta f_{17} =
0$).  The estimator in Eq.~\ref{equation:simple-entropy-estimator} had a standard deviation of $0.070 \pm 0.001$ and a bias of $0.000 \pm 0.002$.
(Uncertainty estimates over replicates were obtained from expressions in Sections II.A and II.B of \cite{ahn-fessler:2003:standard-errors}.)

Recently, an additional estimator for the entropy was
presented~\cite{pappu:jpcb:2010:entropy-estimator} based on
temperature derivatives of the expression for the free energy, and in
the appendix we perform more complete analyscomparison of the esimator
The difference is statistically indistinguishable and suggests that
the simple estimator of Eq.~\ref{equation:simple-entropy-estimator} is
as efficient as the more complex estimator
(Eq.~\ref{equation:rohit-entropy}), as would be expected from the fact
that the MBAR estimates for $f_i$ and $u_i$ are asymptotically
efficient~\cite{shirts-chodera:jcp:2008:mbar}.


%Validation
%* Harmonic oscillators [tests/harmonic-oscillators] JDC
%   Show that estimated free energies and expectations are consistent; uncertainties are correct.  One state has no samples.
%   - relative free energies and entropy/enthalpy decomposition
%   - expectation values of position
%   - expectation values of mean-square displacement (and heat capacity?) show treatment of correlated estimates
%   FIGURE: Confidence interval validation plots for free energies, enthalpies, entropies, position, and mean-square displacement
%   Q: do we do this as a function of number of states to show where it breaks down?
%   A: We can put most of the plots in the Supplementary Info, and only report on, say, the 68% confidence interval in a table in the main paper.  We can try for various numbers of samples, showing at what number of samples the confidence levels are good. 
  
\subsection{Advantages of MBAR over pairwise MBAR}

%Examples
%* Harmonic oscillators with variable spacing [examples/harmonic-oscillators] MRS
%   - illustrate efficiency of sequential BAR vs MBAR in estimating free energy difference between 1 and N
%      + MBAR gives better uncertainties than sequential BAR because of treatment of data reuse     
%   FIGURE: Harmonic oscillators diagram; RMS error vs oscillator spacing for BAR and MBAR; also show average reported uncertainty for BAR and MBAR.

MBAR is guaranteed to always have lower variance than BAR.  However,
theoretical proofs of optimality are not 

In collections of thermodynamic states where there is relatively states have relatively
high overlap, then MBAR is guaranteed to yield free energy differences
with lower variance that BAR.

For example, we choose relatively large overlap with $k=0.2$ in all
cases, and look at the variance of the free energy beween

Of course, for optimally chosen states.  However, in this case,
sequential application of BAR on neighboring states does not correctly
predict the variance, as BAR will include samples from

In some cases, the relationship between neighboring thermodynamics
states might be more complicated than simply a linear chain where
states $i-1$ and $i$ are neighbors, $i$ and $i+1$ are neighbors, and
so forth.  For example, if we perform a series of umbrella sampling
simulations with restraint centers in an $(x,y)$ grid, then the point
$(x_i,y_i)$ has four nearest neighbors $(x_{i-1},y_i)$,
$(x_{i+1},y_i)$, $(x_i,y_{i-1})$ , and $(x_i,y_{i+1})$.  Since the
free energy is a state variable, it is not clear which pairwise path
can be used connect any two states $(x_i,y_i)$ and $(x_j,y_j)$ with
minimum variance.  With MBAR, all states are implicitly considered
simultaneously in a minimum variance way.

Other examples of 2D state with nonlinear connection topologies
include the data collected from simultaneous combination of
Hamiltonian and temperature replica, or Hamiltonian exchange with
multiple coupling parameters.  In all of these cases, it is not
obvious how to best apply pairwise BAR.  



