Monte Carlo统计误差估计

整理自M. E. J. Newman & G. T. Barkema Monte Carlo Methods in Statistical Physics

本文所有时间单位为MCS (Monte Carlo steps per particle),或称为sweep。以single flip为例,进行particle number次的single flip attempt在统计上平均每个particle被尝试翻转一次,定义为一个MCS。其他算法得到的时间需要进行换算rescale。

基本原则

对于$n$个独立的测量$m_i,\ i=1,2,\ldots,n$,标准差为

(独立指测量的MCS间隔为两倍自关联时间以上。)

对于MCMC获得的存在时间关联的测量(意味着不独立),误差与自关联时间$\tau$有关,

这里$t_\text{max}$指除去equilibration后的MCS。

用这样的方法可以估计一次观测量如$E$、$M$的误差。

自关联时间$\tau$的估计可以通过平衡后观测量的自关联函数获得,这里不详细说明。因为是误差的估计,所以自关联时间不需要计算得特别精确。

热容$C$、磁化率$\chi$等一次观测量衍生出的二次观测量,本身需要一组独立测量的一次观测量来计算。以上方法不适用。

Blocking法

Blocking法一种通用的误差分析法。以热容的计算和误差分析为例:

  • 通过MCMC获得了一段$t_\text{max}$长的平衡后的$E$,其自关联时间估计为$\tau$;

  • 将$t_\text{max}$个$E$平分为$n_b$段($n_b$个block),要求每一段的长度要远大于自关联时间,从而能够将每一段看做单独的一次实验。在每一段上各自计算能量期望和热容,然后利用式$\eqref{eqn:1}$(其中$n$换为$n_b$)求标准差。一般取$n_b=10$。如果每一段长度都远大于$\tau$,那么$n_b$的取值应该不会影响误差的数量级。

这个方法简单但粗糙。下面两个方法更可靠一些。

Bootstrap法

一种重采样方法。仍然以热容为例:

  • 从MCMC采样获得的$t_\text{max}$个$E$中可放回地随机抽取$t_\text{max}$个。这$t_\text{max}$原则上应当是独立的,但是略微不独立对boostrap影响也不大。

  • 以这个重采的样本计算热容。

  • 重复以上步骤$n_b$次,得到$n_b$个热容。

  • 计算热容的标准差

    特别注意这里没有式$\eqref{eqn:1}$中的$1/(n-1)$。(一个角度是显然过多的重采样不会减少真的误差。)

虽然bootstrap方法可以容忍一定程度的自关联,但是对于太显著的自关联时间也是无法接受的。个人感觉至少要$\tau/t_\text{max}<0.001$。

一般重采样次数100到1000次。

附pythonbootstrap代码

bootstrap
1
2
3
4
5
6
7
8
9
import numpy as np
from astropy.stats import bootstrap
from astropy.utils import NumpyRNGContext

def boot(x):
n_resample=1000
with NumpyRNGContext(1):
bootresult = bootstrap(x,n_resample,bootfunc=np.mean,)
return (np.mean(bootresult),np.std(bootresult))

Jackknife法

仍然以计算热容为例:

  • 首先从全部$t_\text{max}$中隔$2\tau$选取$E$值,获得长度为$n$的样本,样本中每个测量都是独立的;

  • 从独立样本计算一个热容$C$,作为热容的计算值;

  • 从独立样本中移除第i个$E$,剩下$n-1$个$E$,以此为新的样本计算热容$C_i$;

  • 重复上一步$i=1,2,\ldots,n$(每次得到一个$n-1$长度的样本),共获得$n$个$C_i$。

  • 计算标准差

如果$n>100$,bootstrap的计算量更小,否则jackknife更合适。