Matlab实现大数定理的应用

3年前浏览2195

过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，群号：927550334

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，明天学一点正态分布的评估方法，菲尔兹将指日可待。