2022年数学实验-：数据的统计与分析 .pdf

《2022年数学实验-：数据的统计与分析 .pdf》由会员分享，可在线阅读，更多相关《2022年数学实验-：数据的统计与分析 .pdf（10页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、实验 10：数据的统计与分析习题 5：炮弹射击的目标为一圆形区域，半径为100m,弹着点以圆心为中心成二位正态分布，设在密度函数式当中，x=80m,y=50m,相关系数r=0.4,求炮弹命中圆形区域的概率。1 模型建立设目标中心为坐标原点。Rad(radium)=100,则圆形区域可以表示为：222:xyrad着弹点符合二维正态分布，记其坐标为（x,y）,其概率密度为有：22222211221211(,)exp(2)2(1)21xxyyp x yrrr（1）其中1=x，2=y，由于中心在原点，所以上式中不含有期望值（=0）。于是炮弹命中圆形区域的概率可以利用二重积分求得：(,)Pp x y d

2、xdy（2）以上积分无法用解析访法求解，可以根据Monte Carlo 方法通过下式进行运算：21(2*)(,)(,)mkkkradPp x y dxdyp xyn（3）其中,2(2*)rad表示与圆域外切的正方形区域的面积,n 为投点次数,(,)kkxy表示落在区域中的点的坐标。2 程序设计（程序部分可直接粘贴运行）：1)构造概率密度函数,符合（1）式function f=prob(s1,s2,r,x,y)f=1/(2*pi*s1*s2*sqrt(1-r2)*exp(-1/(1-r2)/2*(x2/s12-2*r*x*y/s1/s2+y2/s22);2)主函数clear alls1=80;s

3、2=50;%s1,s2为标准差r=0.4;n=100000;rad=100;x=unifrnd(-rad,rad,1,n);%在（-100，100）内随机均匀取n组 x,y 值，y=unifrnd(-rad,rad,1,n);sum=0;m=0;名师资料总结-精品资料欢迎下载-名师精心整理-第 1 页，共 10 页 -tic%计时for k=1:n if x(1,k)2+y(1,k)2=rad2%实现 Monte Carlo方法 sum=sum+prob(s1,s2,r,x(1,k),y(1,k);m=m+1;%sum为（3）式右端和式部分 end end toc p=(2*rad)2/n*s

4、um%根据（3）式计算概率3 运行结果及分析：n=100000 1 2 3 4 5 计算结果0.6962 0.6977 0.6965 0.6980 0.6967 计算时间（s）2.161175 2.160833 2.196580 2.142944 2.189436 n=10000 1 2 3 4 5 计算结果0.6906 0.6997 0.6953 0.6972 0.7025 计算时间（s）0.21761 0.20825 0.20871 0.24071 0.22048 最终结果为0.7 左右。通过上表还可以看出，随即试验的次数并不能完全的决定最终结果的准确性。当n=1e5时，其结果比起n=1e

5、4 的结果相对稳定，但是计算时间是后者的10 倍，可以推断若将本方法应用于更大规模的数据处理当中，必然产生精度和计算速度的矛盾。以上问题是实际上反映了局部抽样中必然存在的问题，Monte Carlo 算法的理论基础是Bernoull 大数定理，即：n 次独立重复试验中A 发生的次数k,与 A 在每次试验中发生的概率p 有如下关系：lim|1nkPpn（）而实际中的试验次数必然是有限的，所以最终得到的结果必然会不能完全符合概率值。但是，在多次重复试验中，同分布的随机变量，其总体期望和方差为：_11niiE xExEXn（）_211niiDXD xDxnn（）可以看出，随着试验次数的增加，总体期望

6、并没有发生变化，但是方差变小了。这也就是n=10000 时得到的结果波动性比n=100000 时要强的原因了，试验次数越多，试验结果偏离实际概率的程度就越小。可知，在更大的情况下，对应着一个精度，在该精度要求下，最终结果可以认为是和概率值完全符合。这里不再继续进行次数更多的实验。4 一个错误的分析：同课本中例不同的是，本题的、是相互关联的，即满足二维正态分布，而例种的、坐标是相互独立的，各自满足一维正态分布。因此，在平面上，对两个坐标的取点就必须考虑其相互影响，、的取值不是关于坐标轴对称的（实际是关于原点对称的），因此在计算积分时区域位于四个象限内的积分也不是完全相等的。若此时仍采用例名师资料

7、总结-精品资料欢迎下载-名师精心整理-第 2 页，共 10 页 -算法：用第一象限1的积分值作为整体积分值，就会产生错误（得到的结果为）；类似的，若用第二象限2的积分值作为整体积分值，也会产生错误（得到的结果为）。这都是忽略了、的相互作用造成的。相关系数的意义是，当=时，、完全不相关；r=1 时，、成线性关系；r=1 时，、成负线性关系。用以下程序，绘制不同值时的、prob（x,y）三维图像：r=0.9;%r=0.1/-0.5;x=-100:0.1:100;y=x;X,Y=meshgrid(x,y);Z=prob(s1,s2,r,X,Y);mesh(X,Y,Z)r=0.9 名师资料总结-精品资

8、料欢迎下载-名师精心整理-第 3 页，共 10 页 -r=0.1 r=-0.5 可以看到,三维图形也基本符合钟形分布，但是随着值的变化，概率密度峰值的出现范围随之改变，红色椭圆区域的长轴方向，基本上与走向一致。时，可以推断，概率密度值讲关于、轴对称，即不相关情况（例），可以用某一个象限的值计算。5 另一个方法：利用 MATLAB提供的生成符合二维正态分布的随机二维向量的函数，可以直接生成n个符合题目要求分布的点的x、y 坐标，直接计算生成的点位于以圆点为圆心，100m 为半径的圆内的频率P,根据大数定理（4）式可以得知，当n 趋近于无穷时，频率P就是概率。名师资料总结-精品资料欢迎下载-名师精

9、心整理-第 4 页，共 10 页 -s1=80;s2=50;%初始条件若干，同前r=0.4;n=100000;rad=100;m=0;mu=0,0;%期望sigma=s12,s1*s2*r;s1*s2*r,s22;%协方差矩阵x=mvnrnd(mu,sigma,n);%生成服从二维分布的随机二维向量for k=1:n%检验在圆域内的点的个数if x(k,1)2+x(k,2)2=rad2 m=m+1;endendP=m/n 结果：P=0.69535000000 0.69969000000 0.69737000000 0.69934000000。可以看到最终的结果仍然在0.7 左右。采用本方法，实

10、际上是完全模拟了现实的投弹过程，是一种直接符合大数定理形势的Monte Carlo 方法。习题：轧钢有两道工序：粗轧和精轧，粗轧钢坯时由于各种随机因素的影响，得到的钢材长度成正态分布，其均值可由轧机调整，而方差是设备精度决定的，不能改变；精轧时将轧得到钢材轧成规定的长度（可以认为没有误差）。如果粗轧后的钢材长度达与规定长度，精轧时要把多的部分轧掉，造成浪费；如果粗轧后的钢材长度已经小雨规定长度，则整根报废，浪费更严重。问题是已知的钢材规定的长度和粗轧后的钢材长度的均方差，求可以调整的粗轧时刚才长度的均值，失踪的浪费最小。从以下两种目标函数种选择一个，在l=2m,=20cm 条件下求均值：（）每

11、粗轧一根刚才的浪费最小（）没得到一根规定长度的钢材浪费最小模型建立本题需要建立反映钢材浪费程度的目标函数，并使其最小，但由于涉及到正态分布的概率密度函数，是一个非线优化问题。之后的建模并没有严格按照优化问题的步骤进行，而是采用了更简便的数值扫描的办法直接找到最小点。I.浪费程度可以直接用浪费的钢材的长度表征，建立关于浪费长度的目标函数：由于粗轧长度的不同会造成两种浪费模式，因此对两种模式分别进行研究：1）粗轧的得到的钢材长度小于规定长度，全部浪费；2）粗轧的得到的钢材长度大于规定长度，大于规定长度的部分被浪费；对于模式1，浪费量的“期望”值（实际是不同浪费长度用其概率密度加权后的和）可用以名师

12、资料总结-精品资料欢迎下载-名师精心整理-第 5 页，共 10 页 -下方法求得：1()lwxp x dx(1)x 为粗轧得到的钢材长度，w1 表示模式 1 的浪费总长度的估计值，p(x)粗轧的概率密度函数；对于浪费模式2，由于是部分浪费，所以浪费长度由x 本身变为x-L，且积分区域也将变为L 右方，有下式：2()()lwxl p x dx(2)本题当中，x 服从正态分布，即:221()()exp22xmp x(3)其中 m 为本题所求的正态分布期望m，为方差。写成优化问题的一本形式：22min12.1()2()()1()()exp223llwwst wxp x dxwxl p x dxxmp

13、 xlml（4）最后的不等式约束条件的原因是：m 为正态分布的期望，若其小于标准长度L,很明显将至少有50%的钢材由于粗轧后小于 L 而被直接浪费，这显然不是最优的方法，所以有Lm;考虑到正态分布的3法则，可以认为位于m 点左方距离大于3的点，其概率密度极小，分布函数值（概率）接近于0，若 L 位于该区域，则浪费模式1 出现的概率基本为0，失去了讨论的价值，因此有 mL+3。w1+w2 的形势可以利用积分性质进行化简：名师资料总结-精品资料欢迎下载-名师精心整理-第 6 页，共 10 页 -12()()()()()*1()*1()lllwwxp x dxxl p x dxxp x dxlp x

14、 dxElF lmlF l（5）其中 E 为正态分布的期望，就是m 的值，F(L)表示在 x=L 点的分布函数值。将（5）代入（4），就可以直接采用扫描m 的办法找到w1+w2 的最小值点。II.第二问是一个条件期望的问题，在第一问的基础上，利用条件期望公式可以直接得到：()(|)()E ABE A BP B（6）其中(|)E A B表示在事件B 发生的条件下A 的期望。对于本题，（6）式的意义是：()()()EEP每粗轧一根钢材的浪费每得到一根规定长度钢材的浪费每得到一根规定长度的钢材因此，可以直接利用第一问的结果除以每得到一根规定长度的钢材的概率即可。2程序设计1）第一问clear v=1

15、;for q=2:0.001:2.6%在lml+3*sigma区间内扫描 m值 m=q;s=0.2;l=2;Fl=normcdf(l,m,s);%求F(l)p(v,:)=m,m-l*(1-Fl);%(5)式 v=v+1;%用p(v,:)记录结果,输出 m,w1+w2 endp plot(p(:,1),p(:,2)%绘制浪费期望值p 同粗扎期望值m的关系曲线2）第二问：clear v=1;名师资料总结-精品资料欢迎下载-名师精心整理-第 7 页，共 10 页 -for q=2:0.01:2.6;m=q;s=0.2;l=2;Fl=normcdf(l,m,s);p(v,:)=m,(m-l*(1-Fl

16、)/(1-Fl);%式（6），除以每得到一根规定长度钢管的概率，即：L点（规定长度）的分布函数值v=v+1;endp 3运行结果1）第一问m W1+w2 m W1+w2 2 1 2.335 0.4289 2.001 0.997 2.336 0.429 2.002 0.994 2.337 0.429 2.003 0.991 2.338 0.429 2.004 0.988 2.339 0.4291 2.34 0.4291 2.321 0.4295 2.341 0.4292 2.322 0.4294 2.342 0.4293 2.323 0.4293 2.343 0.4293 2.324 0.429

17、2 2.344 0.4294 2.325 0.4292 2.345 0.4295 2.326 0.4291 2.346 0.4296 2.327 0.429 2.347 0.4297 2.328 0.429 2.348 0.4299 2.329 0.429 2.349 0.43 2.33 0.4289 2.331 0.4289 2.597 0.5998 2.332 0.4289 2.598 0.6008 2.333 0.4289 2.599 0.6017 2.334 0.4289 2.6 0.6027 表 1：数据扫描结果名师资料总结-精品资料欢迎下载-名师精心整理-第 8 页，共 10 页

18、-22.12.22.32.42.52.62.72.82.930.40.50.60.70.80.911.1mw1+w2m=2.333图 1：扫描趋势曲线可以看到，最终的结果m 取 2.33 左右时，可以使得浪费的期望值w1+w2 最小(0.4389m).由曲线可以直观的看到变化趋势，结合本题的实际意义很好理解。以m=2.33 作为起始点，当 m 减小小时，L(始终位于 m 左侧，见”模型建立”)靠近 m,则产生第一类浪费模式（整根报废）的钢材量增加；若m 增大，L 远离 m,第一浪费模式的钢材减少，但第二类浪费模式（精轧浪费）的数量增加。由于精轧过程的浪费一定会少于整根报废的情况，所以曲线右端部

19、分相对于左半部分比较平缓。2）第二问m m 2 2 2.356 0.4479 2.001 1.9861 2.357 0.4479 2.002 1.9723 2.358 0.4479 2.003 1.9586 2.359 0.4479 2.004 1.9451 2.36 0.448 2.361 0.448 2.348 0.4482 2.362 0.448 2.349 0.4481 2.363 0.4481 2.35 0.4481 2.364 0.4482 2.351 0.448 2.352 0.448 2.597 0.6007 2.353 0.4479 2.598 0.6016 2.354 0.4479 2.599 0.6026 2.355 0.4479 2.6 0.6035 名师资料总结-精品资料欢迎下载-名师精心整理-第 9 页，共 10 页 -表 2：数据扫描结果2 最终结果 m 在 2.352.36 之间。可以看出不同的最优条件对结果会产生不可忽略的影响，实际问题中应该根据自身情况制定最合理的最优条件，以次得到合理的决策。名师资料总结-精品资料欢迎下载-名师精心整理-第 10 页，共 10 页 -