Ferrenberg, Alan M., and D. P. Landau. “Critical behavior of the three-dimensional Ising model: A high-resolution Monte Carlo study.” Physical Review B 44.10 (1991): 5081.

在MC计算Ising这样的模型的临界行为时候，需要算法精度非常高。通常计算相变点需要使用有限尺度标度(finite size scaling, FSC)方法，通过求高阶矩对温度曲线关于尺度的不动点来求相变温度。然而，普通的MC算法（SF和WC等）受计算资源限制，不可能无穷地增加计算步数来提高精度，所以几乎不可能达到FSC的精度要求。

Histogram方法是一种减少总体计算量的方法，只需要得到一个温度$T$下的足够大的系综，就可以使用histogram reweighing得到邻近的其他$T$。

原理

单直方图(SH)方法

一般来说，MC在$T=T_0$处产生系统的一系列构型（轨迹），每个构型出现的频率正比于Boltzmann权重$\exp[-\beta_0 H]$，这里$\beta_0=1/k_BT$，$H$是系统的哈密顿量。一种惯用的标记是记$K=-J/k_BT$，成为耦合强度，下文我们主要使用这个记号。这样Boltzmann权重记为$\exp[K_0 E]$，$E$是构型的能量。那么在$K_0$时观测到能量$E$且磁矩$M=\sum_i\sigma_i$的构型的概率为

$\begin{equation}\label{eqn:p0} P_{K_0}(E,M)=\frac1{Z(K_0)}W(E,M)\exp[K_0 E] \end{equation}$

这里$W(E,M)$是有着能量为$E$、磁矩为$M$的构型的总数（简并度、态密度），$Z(K_0)$是配分函数。因为MC模拟都是按照这个概率来产生构型，所以相当于产生一个直方图$H(E,M)$以估计平衡时的概率分布。$H(E,M)$表示采样时能量为E磁矩为M的构型数目。当MC步数无穷多时，频率分布就是概率分布，但是对于有限的MC步数，采样不足就会造成统计误差。然而，尽管误差存在，$H(E,M)/N$仍然是一个对$P_{K_0}$的估计（N是MC的步数）。按照这个思路，我们类比式$\eqref{eqn:p0}$写出

$\begin{equation}\label{eqn:h} H(E,M)=\frac{N}{Z(K_0)}\widetilde{W}(E,M)\exp[K_0E] \end{equation}$

这里$\widetilde{W}(E,M)$是对真实态密度$W(E,M)$的统计估计。

下面我们要证明，如果知道了某一个$K_0$下的概率分布，那么我们可以确定任何$K$下的概率分布。对一般$K$，概率分布为

$\begin{equation} P_{K}(E,M)=\frac1{Z(K)}W(E,M)\exp[K E] \end{equation}$

因为$K_0$下的MC给出了$H(E,M)$，所以这个式子给出了$\widetilde{W}(E,M)$

$\begin{equation} \widetilde{W}(E,M)=\frac{Z(K_0)}{N} H(E,M)\exp[-K_0E] \end{equation}$

用$\widetilde{W}(E,M)$估计$W(E,M)$，那么

$\begin{equation}\label{eqn:mh1} P_k(E,M)=\frac{H(E,M)\exp[\Delta KE]}{\displaystyle{\sum_{E,M} H(E,M)\exp[\Delta KE]}} \end{equation}$

这里$\Delta K=K-K_0$。之后，从$P_K(E,M)$出发可以得到关于$E$和$M$的任何函数的期望值

$\begin{equation}\label{eqn:mh2} \la f(E,M)\ra_{K}=\sum_{E,M}f(E,M)P_K(E,M) \end{equation}$

式$\eqref{eqn:mh1}$和$\eqref{eqn:mh2}$称为单直方图方程（SH方程）。

由于$K$可以连续变化，直方图方法非常适用于确定热力学函数导数的峰。我们可以使用式$\eqref{eqn:mh2}$来计算热力学量的期望以及关于$K$的任意阶导数。例如$\la|m|\ra$关于$K$的一阶导数是

$\begin{equation} \frac{\pd\la|m|\ra}{\pd K}=\la|m|E\ra-\la|m|\ra\la E\ra, \end{equation}$

二阶导数是

$\begin{equation} \frac{\pd^2\la|m|\ra_K}{\pd K^2}=\left(\la|m|E^2\ra_K-\la|m|\ra_K\la E^2\ra_K\right)-2\la E\ra_K\frac{\pd\la|m|\ra_K}{\pd K} \end{equation}$

当一个热力学量在某$K$处达到峰值时，其导数为零。利用这一点，对导函数数值求根，就可以精确计算相变温度。

局限性

因为$H_0$毕竟存在统计误差，实际上能够通过SH方法得到的、远离$K_0$的$K$下的概率分布会逐渐显示出误差。超出了某个$\Delta K$范围，概率分布就不可信，这个$\Delta K$会随着系统尺度增加而减小。

需要注意的是，histogram方法本身不能提高单个温度下物理量的精确程度。初始的温度下仍然需要大量的计算。可以粗糙地想象，初始温度下的计算越精确，那么histogram可用的温度领域就越大。