\section{Alternate Entropy Estimators}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BODY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Recently, Pappu and coworkers proposed a new estimator for entropy differences between thermodynamic states~\cite{pappu:jpcb:2010:entropy-estimator} based on the multistate Bennett acceptance ratio~\cite{shirts-chodera:jcp:2008:mbar}.
The identity that forms the basis of their proposed estimator, given in Eq.~33 of \cite{pappu:jpcb:2010:entropy-estimator}, is written in reduced (dimensionless) units as
\begin{eqnarray}
\lefteqn{s_i = - f_i + \sum_{j=1}^K N_j \expect{\mu_i v_i}_j} \nonumber \\
&& \mbox{} + \left[ \sum_{j=1}^K N_j \frac{T_i}{T_j} (\expect{\mu_i v_j}_j - \expect{\mu_i}_j \expect{v_j}_j)  \right. \nonumber \\
&& \left. + \sum_{j=1}^K \sum_{k=1}^K N_j N_k \frac{T_i}{T_k} \left[ (f_k + s_k) \expect{\mu_i \mu_k}_j - \expect{\mu_i \mu_k v_k}_j \right] \right] \nonumber
\end{eqnarray}
where $s_i$ denotes the entropic contribution to the free energy in dimensionless units (having subsumed the temperature $T_i$ written in Eq.~33 of \cite{pappu:jpcb:2010:entropy-estimator}), $f_i$ is the dimensionless free energy, $N_i$ denotes the number of samples from state $i$, $v_i(x)$ is the reduced potential energy of state $i$, and $\mu_i(x)$ is the weighting function
\begin{eqnarray}
\mu_i(x) &=& \frac{\exp[f_i - v_i(x)]}{\sum\limits_{k=1}^K \exp[f_k - v_k(x)]} \nonumber
\end{eqnarray}
Replacing expectations $\expect{A}_j$ with empirical estimators $N_j^{-1} \sum_{n=1}^{N_j} A(x_{jn})$, we obtain the estimator for $s_i$:
\begin{eqnarray}
\lefteqn{\hat{s}_i = - \hat{f}_i + \sum_{j=1}^K \sum_{n=1}^{N_j} \mu_i(x_{jn}) v_i(x_{jn})} \label{equation:rohit-entropy} \\
&& \mbox{} + \left[ \sum_{j=1}^K  \frac{T_i}{T_j} \sum_{n=1}^{N_j} \mu_i(x_{jn}) v_j(x_{jn}) \right. \nonumber \\
&& \mbox{} - \sum_{j=1}^K \frac{T_i}{T_j} N_j^{-1} \left( \sum_{n=1}^{N_j} \mu_i(x_{jn}) \right) \left( \sum_{n=1}^{N_j} v_j(x_{jn}) \right)  \nonumber \\
&& \mbox{} + \sum_{j=1}^K \sum_{k=1}^K N_k \frac{T_i}{T_k} (\hat{f}_k + \hat{s}_k) \sum_{n=1}^{N_j} \mu_i(x_{jn}) \mu_k(x_{jn}) \nonumber \\
&& \left. - \sum_{j=1}^K \sum_{k=1}^K N_k \frac{T_i}{T_k} \sum_{n=1}^{N_j} \mu_i(x_{jn}) \mu_k(x_{jn}) v_k(x_{jn}) \right] \nonumber
\end{eqnarray}
Because $\hat{s}_i$ appears on both the left-hand side (LHS) and right-hand side (RHS), these equations must be solved self-consistently for the $\hat{s}_i$; this is done by fixing the $\hat{f}_i$ determined from MBAR and iteratively updating the estimate for $\hat{s}_i$ on the LHS~\cite{pappu:jpcb:2010:entropy-estimator}.
Note that only relative entropies, differences of the form $\Delta s_{ij} \equiv s_j - s_i$, are meaningful, just as only relative free energy differences $\Delta f_{ij} \equiv f_j - f_i$ are can be estimated because MBAR provides estimates of $f_i$ only up to an arbitrary additive constant.

The reduced free energy $f_i$, average reduced energy $u_i = \expect{v_i}_i$, and reduced entropy $s_i$ obey the Helmholtz relationship
\begin{eqnarray}
f_i &=& u_i - s_i \nonumber
\end{eqnarray}
An alternative 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}
The free energy $f_i$ (up to an additive constant across all states) can easily be estimated with MBAR using Eq.~11 of \cite{shirts-chodera:jcp:2008:mbar}, and $u_i$ can be estimated using Eq.~15 of the same reference.  
This then provides an alternative, simpler approach to the estimation of the reduced entropy $s_i$ over Eq.~\ref{equation:rohit-entropy}.
In fact, we recognize that the estimator on the RHS of Eq.~\ref{equation:simple-entropy-estimator} corresponds to the first two terms on the RHS of Eq.~\ref{equation:rohit-entropy}, as $\hat{u}_i \equiv \sum_{j=1}^K \sum_{n=1}^{N_j} w_i(x_{jn}) v_i(x_{jn})$.
The remaining terms (those enclosed in brackets), might be surmised to correspond to a quantity that, once the iterative estimation of these equations has converged, is identically or nearly zero.

To test both entropy estimators, we consider a system consisting of a simple one-dimensional harmonic oscillator, for which all thermodynamic properties can be computed exactly.
The system has reduced potential
\begin{eqnarray}
v(x) &=& \frac{\beta K}{2} (x - x_0)^2 \nonumber
\end{eqnarray}
with $\beta = 1$ and $K = 0.2$.
For each experiment, 40 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:rohit-entropy} (Eq.~30 of Ref.~\cite{pappu:jpcb:2010:entropy-estimator} had a standard deviation of $0.071 \pm 0.001$ and bias of $0.001 \pm 0.002$, while 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}.)
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}.

The bracketed terms in Eq.~\ref{equation:rohit-entropy} have a mean of $-0.0004 \pm 0.0003$, suggesting that while the expectation of the quantity in brackets is zero.
While it would be difficult to prove that the bracketed terms increase the asymptotic variance of the estimator over that of the simple estimator in Eq.~\ref{equation:simple-entropy-estimator}, the terms have a standard deviation of $0.0108 \pm 0.0002$, over the replicates due to statistical or numerical noise, suggesting this contributes some small amount to an increase in variance of the estimator in Eq.~\ref{equation:rohit-entropy}.

We also compared the two entropy estimators in the estimation of the entropy change in transferring 3-methylindole from vacuum to water in explicit solvent.
A set of 11 alchemical intermediates were used in which the interactions of the solute and solvent were decoupled in stages, and 50 ns trajectories in the isothermal-isobaric ensemble were generated, with potential energies at all alchemical intermediates stored every 1 ps.
MBAR was applied to the entire dataset to estimate the \emph{target} entropy difference between the vacuum and solvated endstates.
The dataset was then divided into 100 sequential blocks of 500 snapshots each, and both entropy estimators were applied to each block.
The standard deviation and bias of each estimator were computed from the deviation from the \emph{target} entropy difference over the blocks.
Here, the bias of the estimator in Eq.~\ref{equation:rohit-entropy} was significant ($-1.6 \pm 0.4$) compared to that of the simpler estimator of Eq.~\ref{equation:simple-entropy-estimator} ($0.0 \pm 0.3$).  
The standard deviation of the estimator of Eq.~\ref{equation:rohit-entropy} ($-4.0 \pm 0.3$) was also larger than that of the simpler estimator of Eq.~\ref{equation:simple-entropy-estimator} ($3.1 \pm 0.2$). 
%MRS: this doesn't make sense -- you have a negative standard deviation.

Together with the extra complexity of including these terms, this suggests that users of MBAR interested in entropy estimates should prefer the simpler estimator of Eq.~\ref{equation:simple-entropy-estimator} over the more complex estimator of Eq.~\ref{equation:rohit-entropy}.

We also note that the statistical error in entropy differences estimated through the simple expression in Eq.~\ref{equation:simple-entropy-estimator} can be straightforwardly obtained from the covariance estimation equation given in Eq.~10 of \cite{shirts-chodera:jcp:2008:mbar}.
As described in Section IV of \cite{shirts-chodera:jcp:2008:mbar}, we first augment the $N \times K$ weight matrix by another $K$ columns corresponding to the reduced potential $v(x)$ as an observable, such that the corresponding $q_{K+i}(x) = v_i(x) q_i(x)$ for column $K+i$.
The corresponding weight column is given by
\begin{eqnarray}
W_{n,K+i} &=& \hat{c}_{K+i}^{-1} \frac{u_i(x_n) \exp[-u_i(x_n)]}{\sum\limits_{k=1}^K N_k \exp[\hat{f}_k - u_k(x_n)]}
\end{eqnarray}
Making note of the fact that, in terms of the log normalization constants $\theta_i = \ln c_i$ defined in \cite{shirts-chodera:jcp:2008:mbar}, we have
\begin{eqnarray}
\Delta f_{ij} &=& - \ln c_j + \ln c_i = e^{\theta_i - \theta_j} \nonumber \\
\Delta u_{ij} &=& \frac{c_{K+j}}{c_j} - \frac{c_{K+i}}{c_i} = e^{\theta_{K+j} - \theta_j} - e^{\theta_{K+i} - \theta_i} \nonumber \\
\Delta s_{ij} &=& \ln \frac{c_j}{c_i} + \frac{c_{K+j}}{c_j} - \frac{c_{K+i}}{c_i} = \theta_j - \theta_i + e^{\theta_{K+j} - \theta_j} - e^{\theta_{K+i} - \theta_i} \nonumber 
\end{eqnarray}
\begin{eqnarray}
\delta^2 \Delta f_{ij} &=& \hat{\Theta}_{i,i} + \hat{\Theta}_{j,j} - 2 \hat{\Theta}_{i,j} \nonumber \\
\delta^2 \Delta u_{ij} &=& \left[\begin{array}{rrrr} \frac{\partial \Delta u_{ij}}{\partial \theta_i} & \frac{\partial \Delta u_{ij}}{\partial \theta_j} & \frac{\partial \Delta u_{ij}}{\partial \theta_{K+i}} & \frac{\partial \Delta u_{ij}}{\partial \theta_{K+j}}\end{array}\right] \tilde{\bf \Theta} \left[\begin{array}{r} \frac{\partial \Delta u_{ij}}{\partial \theta_i} \\ \frac{\partial \Delta u_{ij}}{\partial \theta_j}\\ \frac{\partial \Delta u_{ij}}{\partial \theta_{K+i}}\\ \frac{\partial \Delta u_{ij}}{\partial \theta_{K+j}}\end{array}\right] \nonumber \\
&=& \left[\begin{array}{rrrr}u_i & -u_j & -u_i & u_j\end{array}\right] \tilde{\bf \Theta} \left[\begin{array}{r}u_i\\-u_j\\-u_i\\u_j\end{array}\right] \nonumber \\
\delta^2 \Delta s_{ij} &=& \left[\begin{array}{rrrr} \frac{\partial \Delta s_{ij}}{\partial \theta_i} & \frac{\partial \Delta s_{ij}}{\partial \theta_j} & \frac{\partial \Delta s_{ij}}{\partial \theta_{K+i}} & \frac{\partial \Delta s_{ij}}{\partial \theta_{K+j}}\end{array}\right] \tilde{\bf \Theta} \left[\begin{array}{r} \frac{\partial \Delta s_{ij}}{\partial \theta_i} \\ \frac{\partial \Delta s_{ij}}{\partial \theta_j}\\ \frac{\partial \Delta s_{ij}}{\partial \theta_{K+i}}\\ \frac{\partial \Delta s_{ij}}{\partial \theta_{K+j}}\end{array}\right] \nonumber \\
&=& \left[\begin{array}{rrrr}(-1+u_i) & (1-u_j) & -u_i & u_j\end{array}\right] \tilde{\bf \Theta} \left[\begin{array}{r}(-1+u_i)\\(1-u_j)\\-u_i\\u_j\end{array}\right] \nonumber 
\end{eqnarray}
where 
\begin{eqnarray}
\tilde{\bf \Theta} &\equiv& \left[\begin{array}{llll}
\hat{\Theta}_{i,i} & \hat{\Theta}_{i,j} & \hat{\Theta}_{i,K+i} & \hat{\Theta}_{i,K+j} \\
\hat{\Theta}_{j,i} & \hat{\Theta}_{j,j} & \hat{\Theta}_{j,K+i} & \hat{\Theta}_{j,K+j} \\
\hat{\Theta}_{K+i,i} & \hat{\Theta}_{K+i,j} & \hat{\Theta}_{K+i,K+i} & \hat{\Theta}_{K+i,K+j} \\
\hat{\Theta}_{K+j,i} & \hat{\Theta}_{K+j,j} & \hat{\Theta}_{K+j,K+i} & \hat{\Theta}_{K+j,K+j}
\end{array}\right] \nonumber
\end{eqnarray}
and the asymptotic covariance matrix estimate $\hat{\bf \Theta}$ is computed according to Eq.~8 of \cite{shirts-chodera:jcp:2008:mbar}.
\end{document}
