蒙特卡洛模拟的简单例子

蒙特卡罗模拟法的理论基础是大数定律。其发展得益于计算能力的提升和计算成本的下降。蒙特卡罗模拟法在很多领域都有广泛的应用。

引子

图片来源：Wikipedia

matrix67的《数学里也能耍流氓》一文，其中提到了旋轮线的面积的证明：

车轮在地上旋转一圈的过程中，车轮圆周上的某一点划过的曲线就叫做“旋轮线”。在数学和物理中，旋轮线都有着非常重要而优美的性质。比如说，一段旋轮线下方的面积恰好是这个圆的面积的三倍。这个结论最早是由伽利略（Galileo Galilei，1564-1642）发现的。不过，在没有微积分的时代，计算曲线下方的面积几乎是一件不可能完成的任务。伽利略是如何求出旋轮线下方的面积的呢？

他的方法简单得实在是出人意料：它在金属板上切出旋轮线的形状，拿到秤上称了称，发现重量正好是对应的圆形金属片的三倍。

“称重”是对蒙特卡罗模拟法最形象的解释。

蒙特卡洛模拟

我们首先把旋轮线放在坐标系中：

旋轮线的方程为：

x=r*(θ-sinθ)

y=r*(1-cosθ)

旋轮线内切于一个矩形，并且将这个矩形分为两部分：旋轮线上方平面和下方平面。

假设点在平面上的分布是均匀的，而且所有点的面积、质量相等。则，旋轮线上方平面和下方平面的面积之比等于质量之比，也等于平面内点的个数之比。

如果从矩形区域内随机取点，随着次数的增加，在旋轮线上方平面内的总数与下方平面内的总数之比将趋向一个固定的数值。

这个数值，就是两个平面的面积之比。只要有了这个比例，计算上下两个平面的面积就简单了。

由于旋轮线是轴对称的，所以我们只研究[0,pi]区间内的面积比。

具体过程如下：

• 首先，在[0,2]内取一个随机数作为y，在[0,pi]内取一个随机数作为x，并且保证x和y是完全独立的。
• 然后，根据旋轮线方程，用y求解θ，再用θ计算x′。则(x′ , y)在旋轮线上。
• 最后，比较x和x′。如果x>x′，则随机取到的点(x , y)在旋轮线下平面内，计数变量counter增加1。否则，(x , y)在旋轮线上平面内。
• 重复以上过程N次，当N足够大时，counter/N的值就是上平面与下平面的面积之比。

用Matlab做100000次随机抽样，counter/N最后收敛至0.7502（这是旋轮线下方平面与矩形面积的比例）。

容易得到，旋轮线下方平面的面积为3.0008pi=0.7502*(2*2pi)。圆面积为pi。结论得证。

对于一些特殊分布的问题，只要调整产生随机数的算法就可以了。

随机数是蒙特卡洛模拟的关键

在TB上，用以下程序进行测试，只要K线根数足够多，Pips最终将收敛至0.75。

但是，在通达信上同样的逻辑，

Pips并不能收敛至0.75。

经过分析，TB中同一根K线上的Seed1和Seed2是两个独立的随机数。

但是在通达信中，同一根K线上的Seed1和Seed2是相同的。也就是说，随机数X和Y并不独立，存在线性关系Y=2/pi*X。

并且，这条直线在旋轮线下方平面内，所以，不管我们计算多少次，条件X2<X总是满足，Pips的值恒等于1。

图中Pips对1的偏离，并不是基于大数定理的收敛，而是计算误差。