数学实验10数据的统计与分析

1、实验 10：数据的统计与分析习题5：炮弹射击的目标为一圆形区域，半径为100m,弹着点以圆心为中心成二位正态分布，设在密度函数式当中，=80m, =50m,相关系数r=0.4,求炮弹命中圆形区域的概率。1 模型建立设目标中心为坐标原点。Rad(radium)=100,则圆形区域可以表示为：着弹点符合二维正态分布，记其坐标为（x,y）,其概率密度为有：（1）其中=，=，由于中心在原点，所以上式中不含有期望值（=0）。于是炮弹命中圆形区域的概率可以利用二重积分求得：（2）以上积分无法用解析访法求解，可以根据Monte Carlo方法通过下式进行运算：（3）其中,表示与圆域外切的正方形区域的面积,n

2、为投点次数, 表示落在区域中的点的坐标。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;s2=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

3、=0;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）式右端和式部分 endendtocp=(2*rad)2/n*sum%根据（3）式计算概率3 运行结果及分析：n=10000012345计算结果0.69620.69770.69650.69800.6967计算时间（s）2.1611752.1608332.1965802.1429442.189436n=1000012345计算结果0.69060.69970.69530.697

4、20.7025计算时间（s）0.217610.208250.208710.240710.22048最终结果为0.7左右。通过上表还可以看出，随即试验的次数并不能完全的决定最终结果的准确性。当n=1e5时，其结果比起n=1e4的结果相对稳定，但是计算时间是后者的10倍，可以推断若将本方法应用于更大规模的数据处理当中，必然产生精度和计算速度的矛盾。以上问题是实际上反映了局部抽样中必然存在的问题，Monte Carlo算法的理论基础是Bernoull大数定理，即：n次独立重复试验中A发生的次数k,与A在每次试验中发生的概率p有如下关系：（）而实际中的试验次数必然是有限的，所以最终得到的结果必然会不能

5、完全符合概率值。但是，在多次重复试验中，同分布的随机变量，其总体期望和方差为：（）（）可以看出，随着试验次数的增加，总体期望并没有发生变化，但是方差变小了。这也就是n=10000时得到的结果波动性比n=100000时要强的原因了，试验次数越多，试验结果偏离实际概率的程度就越小。可知，在更大的情况下，对应着一个精度，在该精度要求下，最终结果可以认为是和概率值完全符合。这里不再继续进行次数更多的实验。 4 一个错误的分析：同课本中例不同的是，本题的、是相互关联的，即满足二维正态分布，而例种的、坐标是相互独立的，各自满足一维正态分布。因此，在平面上，对两个坐标的取点就必须考虑其相互影响，、的取值不是

6、关于坐标轴对称的（实际是关于原点对称的），因此在计算积分时区域位于四个象限内的积分也不是完全相等的。若此时仍采用例算法：用第一象限的积分值作为整体积分值，就会产生错误（得到的结果为）；类似的，若用第二象限的积分值作为整体积分值，也会产生错误（得到的结果为）。这都是忽略了、的相互作用造成的。相关系数的意义是，当=时，、完全不相关；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,

7、Y,Z) r=0.9r=0.1r=-0.5可以看到,三维图形也基本符合钟形分布，但是随着值的变化，概率密度峰值的出现范围随之改变，红色椭圆区域的长轴方向，基本上与走向一致。时，可以推断，概率密度值讲关于、轴对称，即不相关情况（例），可以用某一个象限的值计算。5 另一个方法：利用MATLAB提供的生成符合二维正态分布的随机二维向量的函数，可以直接生成n个符合题目要求分布的点的x、y坐标，直接计算生成的点位于以圆点为圆心，100m为半径的圆内的频率P, 根据大数定理（4）式可以得知，当n趋近于无穷时，频率P就是概率。s1=80;s2=50;%初始条件若干，同前r=0.4;n=100000; rad

8、=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左右。采用本方法，实际上是完全模拟了现实的投弹过程，是一种直接符合大数定理形势的Monte Carlo方法。习题：轧钢有两道工

9、序：粗轧和精轧，粗轧钢坯时由于各种随机因素的影响，得到的钢材长度成正态分布，其均值可由轧机调整，而方差是设备精度决定的，不能改变；精轧时将轧得到钢材轧成规定的长度（可以认为没有误差）。如果粗轧后的钢材长度达与规定长度，精轧时要把多的部分轧掉，造成浪费；如果粗轧后的钢材长度已经小雨规定长度，则整根报废，浪费更严重。问题是已知的钢材规定的长度和粗轧后的钢材长度的均方差，求可以调整的粗轧时刚才长度的均值，失踪的浪费最小。从以下两种目标函数种选择一个，在l=2m, =20cm条件下求均值：（）每粗轧一根刚才的浪费最小（）没得到一根规定长度的钢材浪费最小模型建立本题需要建立反映钢材浪费程度的目标函数，并

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

11、费模式2，由于是部分浪费，所以浪费长度由x本身变为x-L，且积分区域也将变为L右方，有下式：(2)本题当中，x服从正态分布，即:(3)其中m为本题所求的正态分布期望m，为方差。写成优化问题的一本形式：（4）最后的不等式约束条件的原因是：m为正态分布的期望，若其小于标准长度L ,很明显将至少有50%的钢材由于粗轧后小于L而被直接浪费，这显然不是最优的方法，所以有L<m;考虑到正态分布的3法则，可以认为位于m点左方距离大于3的点，其概率密度极小，分布函数值（概率）接近于0，若L位于该区域，则浪费模式1出现的概率基本为0，失去了讨论的价值，因此有m<L+3。w1+w2的形势可以利用积分性

12、质进行化简：（5）其中E为正态分布的期望，就是m的值，F(L)表示在x=L点的分布函数值。将（5）代入（4），就可以直接采用扫描m的办法找到w1+w2的最小值点。II.第二问是一个条件期望的问题，在第一问的基础上，利用条件期望公式可以直接得到：（6）其中表示在事件B发生的条件下A的期望。对于本题，（6）式的意义是：因此，可以直接利用第一问的结果除以每得到一根规定长度的钢材的概率即可。2程序设计1）第一问clearv=1;for q=2:0.001:2.6%在l<m<l+3*sigma区间内扫描m值 m=q;s=0.2;l=2; Fl=normcdf(l,m,s);%求F(l) p(

13、v,:)=m,m-l*(1-Fl);%(5)式 v=v+1;%用p(v,:)记录结果,输出m,w1+w2endpplot(p(:,1),p(:,2)%绘制浪费期望值p同粗扎期望值m的关系曲线2）第二问：clearv=1;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)/(1-Fl);%式（6），除以每得到一根规定长度钢管的概率，即：L点（规定长度）的分布函数值v=v+1;endp3运行结果 1）第一问mW1+w2mW1+w2212.3350.42892.0010.9972.3360.4292.002

14、0.9942.3370.4292.0030.9912.3380.4292.0040.9882.3390.42912.340.42912.3210.42952.3410.42922.3220.42942.3420.42932.3230.42932.3430.42932.3240.42922.3440.42942.3250.42922.3450.42952.3260.42912.3460.42962.3270.4292.3470.42972.3280.4292.3480.42992.3290.4292.3490.432.330.42892.3310.42892.5970.59982.3320.42

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