大数定理的小陷阱

过冷水最近这段时间在做积分的学习工作，之前连续分享了好几期的 蒙特卡洛法应用；你所不知道的Monte Carlo形式；重要性抽样方法实例分享 。求积分的问题会不懂？可是就是在下图求积分过程中翻了车；

y=[0 0 0 0 0 0 0 0 1.71 3.24 4.59 5.76 6.75 7.56 8.19 8.64 8.91 9 8.91 8.64  8.19 7.56 6.75 5.76 4.59 3.24 1.71 0]';
x =[ -8 -7 -6 -5 -4.5 -4.3 -4.1 -4 -3.7 -3.4 -3.1 -2.8 -2.5 -2.2 -1.9 -1.6 -1.3 -1 -0.7 -0.4 -0.1 0.2 0.5 0.8 1.1 1.4 1.7 2 ]';
f= fit( x, y,'smoothingspline' );
figure1 = figure;
axes1 = axes('Parent',figure1,'Position',[0.13 0.11 0.8 0.8]);
hold(axes1,'on');
plot(x,y,'LineWidth',3,'Color',[1 0 0]);
ylabel('$f(x)$','FontSize',24,'Interpreter','latex');
xlabel('$x$','FontSize',24,'Interpreter','latex');
 xlim(axes1,[-10 6]);
 ylim(axes1,[-1 14]);
hold(axes1,'off');
set(axes1,'FontSize',16);
annotation(figure1,'arrow',[0.13 0.97],[0.11 0.11],'LineWidth',2);
annotation(figure1,'arrow',[0.58 0.58],[0.11 0.95],'LineWidth',2);

该函数是一个简单函数用符号法、数值积分都可以获得精确解，实际案例函数很复杂，过冷水在此进行了模型简化，是为了说明问题，不要觉得过冷水在“杀鸡牛刀”。

需要用大数定理求积分

用大数定理求如图所示函数的具体表达式则为(N=1000)

x1=linspace(-8,2,1000);
y1=f(x1);
I1=10*sum(y1)./1000
I1 =
   35.9138

注意我们函数图像其实在（-8~-4）区间为0所以也可以等于：

x2=linspace(-4,2,1000);
y2=f(x2);
I2=6*sum(y2)./1000
I2 =
   35.9806

从理论上分析没有问题，实际操作 it is error!  。这是过冷水犯的第一个错误。问题是时间需求是要从f(x)不为零的地方开始积分，过冷水为了省事做整个区间的积分，就是抱着对零积分不影响结果。

实际问题要稍微做做一下变形，要求分段积分。

当然如果你微积分还会一点点自然而然知道：

这个时候过冷水就沙雕了，用大数定理的或分别在（-8，-1），（-1，2）,取了N（1000）个点，所有就有2N个点，所以应用大数定理就是：

x3=linspace(-8,-1,1000);
x4=linspace(-1,2,1000);
y3=f(x3);y4=f(x4);
I3=10*sum([y3;y4])/(2000)
I3 =
   42.8187

再一次 it is error!到现在为止过冷水都没想明白为什么是错的，有对该问题感兴趣的可以留言讨论。“快就是慢，省事就是费事”简单问题还得从原理入手才可靠。

x3=linspace(-8,-1,1000);
x4=linspace(-1,2,1000);
y3=f(x3);y4=f(x4);
I4=7*sum(y3)/1000+3*sum(y4)/1000
I4 =
   35.9588

经过简单对比自然而然会发现

在早期过冷水做计算的时候一直使用的是I₃,直到出现非常明显的错误的时候，一步一步做核查才知道是最基础的公式用错了，之前怀疑是自己代码什么地方敲错了。花了一晚上去检查错误最终排查完后非常有开心なるほど！

再回过头来讲讲为什么I₁（35.9138）≠I₂（35.9806）由于(-8~-4)函数值为零，在该区间取点分占一部分权重，(-8~-4)取点多了 (-4~2)取点自然就少一点 平均值就小一点。通常小数点后两位的误差我们是可以接受的。好巧不不巧过冷水在自己的问题计算误差放大的不合适，该错误可以用重要性抽样来解决。理论没问题，实际操作的时候就有问题了，还是不要省事严格按要求来。

这两处错误让我对积分和大数定理有了更深刻的了解，在简单的问题都要“勤实践，多思考”。总能从中学到新知识。本期的过冷水的分享就这么多，希望都能对大家有帮助，读者在实践中有遇到相似的问题，都可以留言分享，大家共同学习进步。

往期回顾>>>>>>

sym、sum妙用案例

latex()、ploy2sym()、symsum()的妙用

十万个matlab编程问题征集，欢迎来问

matlab科研绘图模板，拿走不谢【精品干货】

互动专区

在matlab爱好者公众号中回复“QQ”，加入公众号专属Q群(非免费)，与更多matlab爱好者一起交流；回复“原创”，加入原创代码共享Q群，小编原创matlab代码任性领！

如需转载，请在公众号中回复“转载”获取授权，未经授权擅自搬运抄袭的，必将追究其责任！