magnetization, the energy of a spin configuration x was given by E ( x ) = − ∑ <k,l> x$_{k}$x$_{l}$ , a sum over all pairs of neighbors [8]. Periodic boundary conditions were used. Simulations were performed for 26 replicas at the temperatures β − 1 = 0 . 5 , 0 . 6 , ..., 3 . 0. Configurations were equilibrated using 250 iterations in which a cluster move, Woolf’s generalization [9] of the Swendsen-Wang algorithm [10], was followed by attempted replica exchanges be- tween adjacent temperature pairs (alternating between { ( β$_{1}$, β$_{2}$ ) , ( β$_{3}$, β$_{4}$ ) , ..., ( β$_{K}$$_{−}$$_{1}$, β$_{K}$ ) } and { ( β$_{2}$, β$_{3}$ ) , ( β$_{4}$, β$_{5}$ ) , ..., ( β$_{K}$$_{−}$$_{2}$, β$_{K}$$_{−}$$_{1}$ ) } ). Equilibration was followed by 10 5 production steps us- ing Monte Carlo trial moves generated by a transition function T ( x $^{′}$| x ) that flips each spin with probability p$_{flip}$ = 1 /N , where N is the total number of spins. The transition function was the same for every replica. After every 100 steps of production, a replica exchange was attempted between a random replica and its higher-temperature neighbor, energies of each configuration were stored, T ( x $^{′}$| x ) was used to sample configurations from 〈O ( x, x $^{′}$) 〉$_{π}$$_{k}$$_{,T}$ (these are distinct from the Markov chain), and transition probabilities between the dif- ferent configurations were calculated using T ( x $^{′}$| x ) = ( p$_{flip}$ ) || x $^{′}$− x |$^{|}$(1 − p$_{flip}$ ) N −|| x $^{′}$− x |$^{|}$, where || x ′ − x || is the number of spin mismatches [3]. 

For this system, we find that both estimators have similar and relatively accurate ef- ficiencies at lower temperatures (Fig. 1). Closer to and above the critical temperature ( k$_{B}$T = 2 . 269), however, the optimized estimator (Eq. 12) has less variance and bias than the original estimator (Eq. 8) with unit weight. Nonetheless, the high-temperature estimates are substantially worse than those at low temperatures. 

Ultimately, the quality of the optimized expression is subject to the same fundamen- tal limitations as Adib described with the original estimator [3]. For the convergence of averages to not be dominated by rare events, (a) most configurations generated by the tran- sition function should fall in typical regions of phase space for the target distribution and (b) most pairs of configurations from different replicas should have substantial density in T ( x $^{′}$| x ). These requirements hold true in the low-temperature regime, where the spins are relatively ordered and a local transition function suffices, but are not met at higher temper- atures, where configurations have disordered spins. As Adib mentioned, convergence may be improved by fine-tuning the simulation procedure or choosing a more suitable transition function. 

Overcoming this sampling limitation is not the intent of the current work. Rather, the purpose is to improve the analysis of a given data set. While the resulting expression (Eq. 

