权权的《Monte Carlo方法》系列预告出来时我已注意到,但由于近来时间有限而直到今晚才有时间细读了第一篇《积分方法》。总体而言写得非常清楚,希望能够坚持继续下去。由于整个Bayesian统计学的根基就在Monte Carlo计算法,我对这个领域一向很感兴趣。但由于在研究工作中尚没有机会使用Bayesian统计学,因此有关的知识都还属间接经验,能够在本论坛探讨这个话题肯定会受益。
在点评之前先推荐几本参考资料,我相信下面这个书单是相当不错的,可惜本人尚无时间深入钻研:
* 对英文著作尚有心理障碍者可以参考一本出色的中文教科书:冯康先生所著《数值计算方法》的第七章《蒙特卡洛方法》(国防工业出版社,1978);
* 一本可读性极强的英文专著,美国哈佛大学教授Jun Liu所著"Monte Carlo Strategies in Scientific Computing" (Springer 2002);
* 对Monte Carlo方法在Bayesian统计学中的广泛应用有兴趣者可以适当参考Andrew Gelman等人所著的"Bayesian Data Analysis" (Second Edition, 2003)之第三部分。
权权将贴子发在物理论坛的目的显然是强调该方法在物理上的运用,我选择在数学论坛加以点评是更看重其统计学背景,着重点各有不同。
>>蒙特卡洛(Monte Carlo)是摩纳哥公国一个城镇,位于地中海沿岸,以其赌场和豪华酒店而闻名,所以就有了以随机方法应用于数值计算的一类方法,被称为Monte Carlo
有关Monte Carlo方法历史背景的最精确描述来自Jun Liu的专著,他指出一批物理学家在二战期间为估算薛定谔方程的本征值而发明了一种基于统计抽样的数值计算法,其最初想法归功于Ulam。后来Ulam的同事Metropolis将该方法命名为Monte Carlo。1950年代Metropolis和几名统计物理学同事发表了一篇经典论文,提出了Markov Chain Monte Carlo(MCMC)算法。而MCMC法后来是Bayesian统计学能够不断前进的主要动力。
>> I = ∫ f(x)p(x)dx
这里可以强调一下x是个矢量。而这个积分是概率统计中数学期望的基本定义,可以写成E(f(x))。对于初学者而言,不要忘记概率密度函数p(x)的取值是可以大于1的,归一化条件是对累积密度函数而言。
>>上述变换就是Monte Carlo积分的基本精神,因为需要用到随机抽样,必然伴随统计误差。
需要用到随机抽样,其动机是想用数值模拟实验中的频率来直接估计一个概率值,而这个概率值是计算许多复杂高维积分的关键。而数值模拟需要产生一个序列的随机数来保证抽样过程的随机性。
>>因为x_i是按照概率密度p(x)分布的随机变量,f(x_i)也是随机变量
为了论述的清晰,应该说x_i是一个随机矢量,那么f(x_i)就是随机变量(标量)。
>>而中心极限定理告诉我们,一组独立随机变量之和的概率分布是高斯,其方差等于每一项随机变量的方差之和
这里关于“中心极限定理”的表述不够精确,容易引起读者混淆,特将Kai-Lai Chung(钟开莱,我国著名数理统计大师许宝禄先生的弟子)概率论教科书中的定义按我的理解方式用英文转述一下:
[Central Limit Theorem] For mutually independent (or weakly correlated) random variables X_1, X_2, ..., X_n with mean mu and variance sigma^2,
√n ( Xbar - mu) / sigma --> N(0,1) in distribution,
where N(0,1) stands for standard Gaussian distribution. This means that the distribution shape of Xbar is more and more like a Gaussian random variable as n increases.
权权的中文表述中漏说了这组随机变量必须来自同一个总体(population)这个重要条件,而且“是高斯”必须改成“在n不断增大时趋向于高斯分布”。
>>而计算高维积分时,Monte Carlo 方法是较优的选择。
权权只从收敛速度的视角来说明Monte Carlo方法在高维情形下的优越性是不够的,更关键的一点是---Monte Carlo模拟结果的精度和概型的维数D无关!结果的精度显然比收敛速度更为重要,因此Monte Carlo方法特别适合求解高维问题。
另外要指出Monte Carlo方法以O(1/√N)的速度收敛,这在理论上已经无法改善。关键要在实际应用中通过巧妙设计模拟概型和改进抽样方法来降低方差。降低方差的技巧是衡量各种Monte Carlo方法优劣的重要指标。
>>不妨回顾一下布丰投针实验来结束本篇,在分布着等距平行木纹的地板上投针,要求针的长度小于木纹之间的距离,几何概形计算结果表明,针与一条木纹相交的概率可以用针的长度、木纹间距和圆周率π表示。而用几何概形计算概率实质上归结为面积的计算,也就是积分的计算,布丰投针实验可以说是用随机抽样计算积分的始祖。
Buffon投针实验看似简单,其中蕴含的几何概型思想值得细细品味。令针的长度为L,木纹间距为S,要求L < S。若针的中点到最近的一条平行线的距离为H,用a表示针与平行线的夹角。显然有约束条件0 <= H <= S/2 和 0 <= a <= π。为了使针与平行线相交,必须满足
H <= (L/2) sin(a)
这样针与平行线相交的概率就是两块面积的比值:
p = ∫_0^π (L/2) sin(a) da / (π S/2 ) = 2L / (π S)
这就是权权所说“而用几何概型计算概率实质上归结为面积的计算,也就是积分的计算”。倘若上式分子中的积分是一个复杂的高维积分,我们就可以用Monte Carlo方法模拟出的p值来估算它。当然假如我们感兴趣的是无理数π值的估算,那么由上式可推出:
π_hat = lim 2L / (S p_n)
极限中的n趋向于正无穷。
希望权权在接下来的系列文章中能谈到以下四种Monte Carlo抽样方法:
* Crude Sampling
* Acceptance-Rejection Sampling
* Stratified Sampling
* Importance Sampling
若能谈及MCMC类方法在统计物理学上的运用则更能引人入胜。