公众号之前有讲了好几期关于Monte Carlo算法的推文。过冷水自以为感觉能够让大家明白什么是Monte Carlo算法。只叹数学方法的深奥灵活岂是一朝一夕就可以掌握的,本期过冷水就和大家分享一下大家所不知道的Monte Carlo算法。
求解定积分:
在被积函数f(x)相当复杂时,就只能采取数值积分的求法。并不是所有的复杂的定积分问题的数值解都可以用决定论算法求解。这个时候就可以用Monte Carlo 方法:在(a,b)区域内均匀随机抽样得到N个点x1、x2、x3、...、xN;求这些点上被积函数的值f(x1)、f(x2)、f(x3)......f(xN)、f(x1)、于是f(x)在[a,b]区域的平均值:
于是定积分:
这数学描述也太简单了吧!求平均数的方法用两个不同公式来表述。两个公式自然而然就联系起来了,用一个公式计算另一个公式的相关参数。这个过程你都感觉不到统计力学的身影,这就是Monte Carlo的另一种思想。平均数的概念都能被玩出花来。数学家群体是一群奇迹般的存在。这个方法和之前所讲的算法完全不一样。
过冷水之前以为的Monte Carlo算法是通过随机撒点求所求区域占规则形状的面积比值然后用规则面积*比值即为所求面积面积。我们来实战演示一下两种思路求积分的具体过程。
Monte Carlo方法一:概率求积分
warning off
feature jit off
n=100000;
x=1+7*rand(1,n);y0=5.*rand(1,n);
y=(x.^4.*exp(x))./(exp(x)-1).^2;
points=find(y0<=y);
k=length(points);
S=7*5*(k/n)
S = 23.2418
Monte Carlo方法二:大数定理求积分
warning off
feature jit off
n=100000;
x=1+7*rand(1,n);y0=5.*rand(1,n);
y=(x.^4.*exp(x))./(exp(x)-1).^2;
I=7*sum(y)/n
I =
23.2552
符号积分法:
warning off
feature jit off
syms x
f=(x.^4.*exp(x))./(exp(x)-1).^2;
y=double((int(f,1,8)))
y =
23.2667
可以看出三种方法计算的积分值相近,这也说明了Monte Carlo的两种不同计算路径得到的计算值都是可信的。
学习的越是深入才能更加清晰的感受到自己所了解的知识的浅薄程度,如果有人问我如何衡量Monte Carlo方法的计算精度,昨天的我会说用平方来比较啊!今天的我就不会这么说,我会给你普及置信区间的概念。和你聊置信水平。这就是日日学习,日日精进。今天学一点Monte Carlo,明天学一点正态分布的评估方法,菲尔兹奖它日可期。