\section{Included code tests}

There are a number of sets of tests included in the pymbar 2.0
distribution at https://simtk.org/pymbar, designed to cover a wide
variety of the uses that users may have for MBAR.

\subsection{harmonic-oscillators}

Tests{\em harmonic-oscillators} test set allows replication of the results in section~\ref{valdiation}, and allows 

for the calculation of free energy differences and expectation values
of position, squared position, mean square discplament,and potential
energy for each oscillator.  Users can set the number of harmonic
oscillators, their equilibrium location, force constant, and number of
samples at from each potential, as well as the overall temperature.
Additionally, users can calculate entropy and energy differences
between the oscillators.  Additionally, the entire experiment can be
conducted multiple time and the error estimates compared with
analytical answers to validate the MBAR formalism.

\subsection{constant-force-optical-trap}

\subsection{parallel-tempering-2dpmf}


\subsection{umbrella-sampling-pmf}

\subsection{alchemical-free-energy}

One common use of MBAR is the calculation of free energies of
alchemical tranformations.  In this test we have included data.  

The test allows the comparison of TI with the trapezoid rule, MBAR,
BAR, and forward and reverse exponential averaging.  The particular
script assumes a free energy format now being used for 


\section{pymbar implememtation}

Further description 

\subsection{helper functions}

Several helper functions are included in the pymbar distribution.
First, the python module {\em timeseries.py} includes a number of free energies

The python module {\em confidenceintervals.py} takes a number of
experiments and compares the distribution of error against an
analytical copy of the variance.

\subsection{solving the equations}

We have implemented both self-consistent iteration and Newton-Raphson
iteration to solve the implicit free energy equations.  Newton-Raphson
can diverge in some cases.  For high precision estimates,
Newton-Raphson is definitely preferred, as self-consistent iteration
can converge extremely slowly near the solution.

The initial guess are default set to zero, but initial guesses can
also be set by performed by calculating faster pairwise BAR estimate
between each pair of states.  Some cases, pairwise overlap may not be
possible to converge even though the full MBAR equations converge
easily -- the umbrella sampling data set included in the pymbar
distriution is an example of this.  In this case, solving directly
from zero is required.


\subsection{expectation values}

%When computing averages, what should we do:

We compute an ensemble average of a property that changes with the
thermodynamic state.  A direct estimate is given by:

\begin{eqnarray}
\expect{A}_i &=& \frac{\int_{\bf{x}} A(x) q_i(x)}{\int_{\bf{x}} q_i(x)} = \frac{C_i}{c_i}
\end{eqnarray}

Where $C_k$ and $c_k$ are the two integrals. However, we can replace
each of these integrals by minimum variance estimate of the integrals
from MBAR, which is:

\begin{eqnarray}
\expect{A}_i &=& \frac{C_i}{c_i} \approx \frac{\sum_{j=1}^K\sum_{n=1}^{n_j} \frac{A(x_n) q_i(x_n)}{\sum_{k=1}^{K} N_k c_k^{-1}q_k{x_n}}}{\sum_{j=1}^K\sum_{n=1}^{n_j} \frac{q_i(x_n)}{{\sum_{k=1}^{K} N_k c_k^{-1}q_k{x_n}}}}
\end{eqnarray}

This can be integrated easily by using $q_A(x) = A_i(x) q_i(x)$ as one
of the states in the standard MBAR implementation, with $N_A = 0$
samples.

For many observables, such as structural observables, $\expect{A}_i$
is a function of state $i$, through the ensemble average weighting,
the actual observable $A(\bfm{x})$ is not; for example, the
instantaneous radius of gyration of a protein is a function of only
the coordinates $\bfm{x}$, not the temperature, though it will
generally increase with temperature as the ensemble of coordinates
$\bfm{x}$ change.

However, if we are sampling from states, then the potential energy is
also a function of the thermodynamic state. If we want the average
value of the potential energy in the $i$th thermodynamic state, then
we need to average:

\begin{eqnarray}
\expect{A}_i &=& \frac{C_i}{c_i} \approx \frac{\sum_{j=1}^K\sum_{n=1}^{n_j} \frac{A(i,x_n) q_i(x_n)}{\sum_{k=1}^{K} N_k c_k^{-1}q_k{x_n}}}{\sum_{j=1}^K\sum_{n=1}^{n_j} \frac{q_i(x_n)}{{\sum_{k=1}^{K} N_k c_k^{-1}q_k{x_n}}}}
\end{eqnarray}

Note that this then requires computing the $i$th thermodynamic state
in the sample from the $j$th state. The pymbar {\bf
  computeExpectations} function supports this functionality by
accepting either a 2D or a 3D array.  In the case of a 2D array, the
first dimension indicates the thermodymic state the observable is
measured from, and the last is the value of the observable at each
sample.  In the case of a 3D array, the first dimension is the state
the observable was sampled from, the second dimension is the label for
the thermodynamic state of the observable, and the last dimension is
the value of the observable at the sample.  

One source of confusion between the two cases can emerge in parallel
tempering, when the observable is the energy.  The energies themselves
do not change with thermodynamic state, but the reduced potentials do.
So if the energies are changed to reduced potentials, a 3D array must
be used; if left as energies, then a 2D array should be used.
However, both cases will yield identical results up to numerical
error.
  
\subsection{Energy and Entropy}

A particularly difficult problem. MBAR can be used to calculate the average
energy at each state, and then subtract this result from the free
energy of that state to obtain the entropy.   MBAR is guaranteed 

An alternate estimator of the entropy was recently presented.
