搜索
数学实验报告实验8数据统计与分析1实验8数据统计与分析分1黄浩2011011743一、实验目的1.掌握概率统计的基本概念及用MATLAB实现的方法。2.练习用这些方法解决实际问题。二、实验内容1.《数学实验》第一版（问题5）问题叙述：与例6类似，但炮弹射击的目标为一半径为100m的圆形区域，弹着点以圆心为中心呈二维正态分布，设在密度函数式中σx=80m,σy=50m，相关系数r=0.4，求炮弹命中圆形区域的概率。建立模型与实验过程：若以原点为圆心，则目标区域：：而炮弹落在任意一点（x,y）的概率密度为：()√{()[]}其中，，下面，我们使用三种方法，分别实现这个计算过程：（一）蒙特卡洛——均值估计法这道题与书中例题十分相似，但不同在于，由于相关系数r的存在，本题不能以第一象限的四分之一圆来替代整个目标区域，事实上，该概率密度函数在一三象限和二四象限分别对称，但在一四象限或者二三象限却不是对称的。于是我直接将整个圆形作为计算区域，因此，炮弹命中圆形区域的概率：∬()(())∑()数学实验报告实验8数据统计与分析2其中，()为落在内的点。那么，我们只要在ϵ[]ϵ[]内随机均匀撒n个点，然后对落在区域内的点进行求和，就能估计出所求的二重积分以及命中概率。所用代码如下：n=1000000;%计算次数sx=80;sy=50;r=0.4;x=unifrnd(-100,100,2,n);%随机撒点k=1/(2*pi*sx*sy*sqrt(1-r^2));z=0;u=0;fori=1:nif(x(1,i)^2+x(2,i)^2)=10000%落在圆内u=k*exp(-0.5/(1-r^2)*(x(1,i)^2/sx^2-2*r*x(1,i)*x(2,i)/sx/sy+x(2,i)^2/sy^2));z=z+u;endendP=40000/n*z重复五次，结果如下：实验序号12345结果0.6981745480.6973596710.697476410.6974708120.699117329均值为：0.6979，标准差为6.65*10^-4由此可见，炮弹命中圆形区域的概率大约为0.7（二）数值积分法在进行完上一方法的实验后，我觉得还可以用数值积分的方法来计算这个二重积分，然而，在之前的实验中，我们仅涉及到了一维的数值积分，没有给出非矩形的二重积分的matlab实现方法。于是我根据积分的最原始意义，将ϵ[]ϵ[]这个区域网格化，每个网格为1*1的正方形，如果网格中心点在内，则代入()，如此不断计算并求和，最后乘以单个网格面积，即为所求概率。所用代码如下：A=ones(200);%对网格矩阵进行预设z=0;sx=80;sy=50;r=0.4;k=1/(2*pi*sx*sy*sqrt(1-r^2));fori=99.5:(-1):-99.5A((100.5-i),:)=-99.5:1:99.5;%矩阵A的元素为相应网格中心点的坐标数学实验报告实验8数据统计与分析3endfori=1:200forj=1:200if(j-100.5)^2+(100.5-i)^2=10000%落在圆内z=z+k*exp(-0.5/(1-r^2)*((j-100.5)^2/sx^2-2*r*(j-100.5)*(100.5-i)/sx/sy+(100.5-i)^2/sy^2));endendendz结果为：0.698061539560864，这与之前的蒙特卡洛算法所得结果基本一致。（三）蒙特卡洛——随机投点法在进行完这步实验后，我又想到，还可以直接设一个二元正态分布X，然后使用rnd命令来生成该满足该分布的随机数，最后统计落在区域内的个数，除以随机数的总数量，便是所求的概率。这种方法的思想是，在重复次数足够多的情况下，以频率代替概率，即书中的“随机投点法”。所用代码如下：n=1000000;sx=80;sy=50;r=0.4;sxy=sqrt(r*sx*sy);mu=[0,0];sigma=[sx^2,sxy^2;sxy^2,sy^2];%期望和协方差矩阵X=mvnrnd(mu,sigma,n);m=0;%定义X为该二元正态分布fori=1:nif(X(i,1)^2+X(i,2)^2=10000)%落在圆内m=m+1;%统计个数endendp=m/n重复五次，结果为：实验序号12345结果0.6976080.6972670.6978610.6981940.697545均值为：0.6977，标准差为：3.13*10^-4得出结论：炮弹命中圆形区域的概率大约为0.7实验感悟：在本次实验中，我使用了蒙特卡洛方法和数值积分两种方法进行计算，而在数学实验报告实验8数据统计与分析4蒙特卡洛方法中，又分为随机投点法和均值估计法两种算法。通过计算和比较，这三种方法的结果都近似为0.7，说明都可以实现概率估计，证明了原理的可行性。但是从算法效率和计算精度上看，这三者是有区别的。数值积分是一种“分割”的思想，从数学思想上说是最原始的，但却是最直观的。在本题中，为了简便，我设置的分块是1*1的正方形，分割较粗，因而精度较低。蒙特卡洛方法是以大数定理为基础。随机投点法，是以频率估计概率，主要应用了伯努利大数定律，而均值估计法，还用到了以平均值估计期望值，即辛钦大数定律。在投点次数相同的情况下，随机投点法所得结果的标准差较小，精度更高，这可能是由于该方法只使用了一个估计，而均值估计法使用了两个估计（频率估计概率、平均值估计期望值）的缘故。但是，均值估计法却巧妙地避开了非均匀的投点过程，可以将概率密度函数的定积分转变为m个概率密度值的求和，在概率分布比较复杂时，由于matlab自身的限制，难以创建一些不常见的分布，因而随机投点法可能完全无法实现，但是均值估计却能顺利进行，这是均值估计法的最大优点。2.《数学实验》第一版（问题7）问题叙述：对于报童问题，如果报纸的需求量服从正态分布𝑁(𝜇，𝜎)，且批发价为𝑎𝐴(𝑛𝐾)，其中𝑛为购进报纸的数量，𝐾为一个给定的常数。建立报童为获得最大利润的数学模型。当已知𝜇，𝜎，𝐴，𝐾，𝑏，𝑐3时，为了获得最大的利润，求解报童每天购进的报纸数量𝑛。建立模型：设每天的报纸需求量为~𝑁(𝜇，𝜎)，报童每天购进报纸𝑛份，购进价为𝑎𝐴(𝑛𝐾)，零售价为b，退回价为c。当𝑛时，报童的利润为(𝑏𝑎)(𝑎𝑐)(𝑛)元；当≥𝑛时，利润为𝑛(𝑏𝑎)元。而r是呈正态分布的，因此将两种情况下将利润与需求概率𝑓()相数学实验报告实验8数据统计与分析5乘并求和即得报童每天的平均利润𝑉(𝑛)，即：𝑉(𝑛)∑[(𝑏𝑎)(𝑎𝑐)(𝑛)]𝑓()𝑛−𝑥=0∑[(𝑏𝑎)𝑛]𝑓()∞𝑥=𝑛此处忽略了x0的情况，这是根据实际问题所做的忽略，当𝜇很大而且𝜎很小时，x0时的概率十分接近于0，可以省略。另一方面，从上式可见，虽然x是离散值，但是由于其值较大，可以将它看做连续随机变量，因而报童每天的平均利润可以写作𝑉(𝑛)∫{(𝑏𝑎)(𝑎𝑐)(𝑛)}()𝑛0∫(𝑏𝑎)𝑛()+∞𝑛∫{(𝑏𝑐)𝐴𝑛(𝑛𝐾)𝑐𝑛}()𝑛0∫{𝑏𝑛𝐴𝑛(𝑛𝐾)}()+∞𝑛其中()√𝜋𝜎𝑒−(𝑥−𝜇)2𝜎2为了求得n使𝑉(𝑛)最大，对𝑉(𝑛)求导，得：𝑉′(𝑛)(𝑏𝑐)𝑛(𝑛)[𝐴(𝑛𝐾)𝑐]∫()𝑛0[𝐴𝑛(𝑛𝐾)𝑐𝑛](𝑛)[𝑏𝐴(𝑛𝐾)]∫()+∞𝑛[𝑏𝑛𝐴𝑛(𝑛𝐾)](𝑛)[𝐴(𝑛𝐾)𝑐]∫()𝑛0[𝑏𝐴(𝑛𝐾)]∫()+∞𝑛令导数为零，可得：∫()𝑛0∫()+∞𝑛𝑏𝐴(𝑛𝐾)𝐴(𝑛𝐾)𝑐当𝜇比标准差𝜎大得多时，可作近似：∫()𝑛0∫()𝑛−∞又因为：∫()+∞𝑛∫()𝑛−∞化简可得：∫()𝑛−∞𝑏𝐴(𝑛𝐾)𝑏𝑐………………(∗)数学实验报告实验8数据统计与分析6而且，当该等式左侧大于右侧时，𝑉′(𝑛)，右侧大于左侧时，𝑉′(𝑛).根据参数的具体情况，对该方程求根（如果存在），即可得最优的报纸购进量。实验过程：根据参数：𝜇，𝜎，𝐴，𝐾，𝑏，𝑐3，先假设(∗)方程有根，而且此根正是V(n)的最大值点。在一开始，我打算使用二分法，即在∈[]的范围内，以5000为初值，不断代入∫()𝑛−∞中并与𝑏−𝐴(−2𝑛𝐾)𝑏−𝑐进行比较，当二者的差值绝对值大于开始规定的最大误差限时，进行二分，直到二者的差值绝对值大于最大误差限、或者恰好分到了零点，停止区间分割，并输出此时的n值。然而，在进行预实验的过程中，发现当n=0时，差值为0，当n为0.01时，差值为负，当n=10000时，差值也为负，这说明在正半轴上，(*)有至少两个根，无法直接使用二分法，不妨先作图来观察该差值的正负变化。使用代码：n1=5;n2=10000;mu=2000;sigma=50;A=0.5;K=50000;b=0.5;c=0.35;n=5000;tol=10^-3;n=0:1:10000;z=normcdf(n,mu,sigma)-(b-A*(1-2*n/K))/(b-c);plot(n,z);绘图：数学实验报告实验8数据统计与分析7（需要注意到，上图中的y轴坐标符号是(*)式的左端减右端所得结果，与𝑉′(𝑛)的符号是相反的。）因此，在区间[0,1969]内，V(x)单调上升，在区间[1969,7455]内，V(x)单调下降，在x7455之后，V(x)单调上升。在n=1969附近有极大值点。但是，当n很大时，虽然可能有比极大值点更高的收益，但每天至少进货7455份，甚至更大，这对一个报童一天的售量来说，是难以承受的。因此，n很大的情况不符合实际，不予考虑。于是我们可以直接在区间[1000,3000]内使用二分法，代码如下：n1=1000;n2=3000;%二分区间mu=2000;sigma=50;A=0.5;K=50000;b=0.5;c=0.35;n=2000;tol=10^-5;z=normcdf(n,mu,sigma)-(b-A*(1-2*n/K))/(b-c);z1=normcdf(n1,mu,sigma)-(b-A*(1-2*n1/K))/(b-c);z2=normcdf(n2,mu,sigma)-(b-A*(1-2*n2/K))/(b-c);while(abs(z)=tol)%差值大于误差限时停止二分if(z*z1)0n2=n;n=(n1+n2)/2;elseif(z*z1)0n1=n;n=(n1+n2)/2;数学实验报告实验8数据统计与分析8else%二分到零点时也停止循环breakendz=normcdf(n,mu,sigma)-(b-A*(1-2*n/K))/(b-c);z1=normcdf(n1,mu,sigma)-(b-A*(1-2*n1/K))/(b-c);z2=normcdf(n2,mu,sigma)-(b-A*(1-2*n2/K))/(b-c);endn结果为：n=1968.2得出结论：每天的进货量为1968份时，可获得最大收益。实验感悟：本题我使用了二分法进行计算，相对于迭代法，这种方法的好处是可以有效减少计算次数，但是，对于此题而言，却完全没有必要用二分法，一是因为计算量很小，二是因为二分法要求端点函数值异号且在该区间内单调，而从得到的图线来看，V(n)在正半轴上是增-减-增，绘图的过程其实就是迭代的一部分了。而且二分法的matlab代码也比较繁琐，因此，以后在处理此类问题时，可以先用小步长搜索合适区间，然后直接使用迭代法即可。三、实验总结本次实验是利用概率论的随机变量分布的知识来进行matlab的代码实现，并解决两个实际问题。第一个问题是求解事件概率，对于连续随机变量而言，即为求解定积分，我使用了蒙特卡洛方法（包括随机投点法和均值估计法）和数值积分法三种方法求解