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

实验10数据的统计与分析化学工程系化71李骥聪2007011861【实验目的】1．掌握概率统计的基本概念及用MATLAB实现的方法；2．练习用这些方法解决实际问题。【实验内容】1.题目3.某校60名学生的一次考试成绩如下：937583939185848277767795948991888683968179977875676968848381756685709484838280787473767086769089716686738094797877635355作直方图，计算均值、标准差、极差、偏度、峰度。解利用MATLAB的计算。输入m文件如下：X=[937583939185848277767795948991...888683968179977875676968848381...756685709484838280787473767086...769089716686738094797877635355];[n,y]=hist(X),%频数表hist(X),%直方图x1=mean(X),%均值x2=median(X),%中位数x3=range(X),%极差x4=std(X),%标准差x5=skewness(X),%偏度x6=kurtosis(X)%峰度解得：n=[203571110958]y=[55.259.664.068.472.877.281.686.090.494.8]x1=80.1000x2=80.5000x3=44x4=9.7106x5=-0.4682x6=3.1529结果：均值为80.1，中位数为80.5，极差为44，标准差为9.7106，偏度为-0.4682，峰度为3.1529。偏度g1反映分布的对称性，g1=-0.46820为左偏态，则样本中，数据位于均值左边的比有变的多。峰度g2反映分布偏离正态分布的程度，g2=3.15293，则样本中，含有较多远离均值的数据，数据不够集中。下图为该学校60名同学成绩的直方图：从最低分53分开始到最高分96分，每隔4.4分为一个区间进行统计，可以看到区间[75,79.4]人数最多为11人，区间[57.4,61.8]人数最少为0。小于均值的人数较多，这与偏度小于0相符合。分数分布不集中，[92.6,96]分的人数较多，即高分段人数偏多，使得分布偏离正态分布。2.题目6.冰淇淋的下部为锥体，上部为半球。设它由锥面22zxy和球面222(1)1xyz围成，用蒙特卡罗方法计算它的体积。模型及其求解因为冰淇淋由锥面和球面222(1)1xyz围成，所以冰淇淋在长方体1,1xy，02z中，若长方体中的点满足222211xyzxy，则该点在冰淇淋内。为了方便编程计算，将冰淇淋z坐标下移一个单位，体积不变。则冰淇淋在长方体1,1xy，11z中，若长方体中的点满足222211xyzxy，则该点在冰淇淋内。长方体的体积为2228，在长方体中任意取点落在冰淇淋内部的概率为p，则冰淇淋的体积为8p用蒙特卡罗方法计算，MATLAB程序如下：n=50000;x=unifrnd(-1,1,3,n);k=0;22zxyfori=1:nif((x(1,i)^2+x(2,i)^2)^0.5-1)=x(3,i)&(x(3,i)^2+x(1,i)^2+x(2,i)^2)=1k=k+1;endendp=k/nV=8*p重复计算4次，计算结果为：p=0.3919p=0.3933p=0.3917p=0.3935V=3.1350V=3.1464V=3.1336V=3.1480当n提高到100000时，重复计算4次，计算结果为：p=0.3933p=0.3939p=0.3927p=0.3930V=3.1467V=3.1511V=3.1412V=3.1440可以看出冰淇淋的体积约为3.14。结果分析下面用微积分方法计算冰淇淋的体积。由二重积分得：2222dd(11)ddddDDDVzxyxyxyxyxy锥面和球面的相交处有222211xyxy，解得：22:1Dxy。所以222222221121212200001123/22300dd(11)ddddd(11)ddd1112[(1)]23231112()233DxyxyVzxyxyxyxyxyrrrrrrrr故冰淇淋体积的精确解为。可以看出用蒙特卡罗方法计算得到的结果具有随机性，精度较低，且在增加取点数目后，精度提高不大。一般来说精度为1/2n阶，当n增加时精度提高较慢。当n=100000时，运算量已经相当大了，所以蒙特卡罗方法用于对非常复杂的定积分、重积分时有较大优势，本题积分较简单，使用蒙特卡罗法意义不太大。本题也是用随机投点法计算的一种办法。3.题目7.对于报童问题，如果报纸的需求量服从正态分布2(,)N，且批发价为(1)naAK，其中为n购进报纸的数量，K为一个给定的常数。建立报童为获得最大利润的数学模型。当已知2000，50，0.5A，50000K，0.5b，0.35c时，为了获得最大的利润，求解报童每天购进的报纸数量n。模型及其求解由题意可知，每天需求量2(,)rN，记报童每天购进报纸的份数为n，当rn时报童售出r份，退回nr份，而每售出1份报童赚ba，退回1份报报童赔ac，所以报童的利润为：2()()()[(1)][(1)]()()()nnnbaracnrbArAcnrbcrAncnKKK当rn时报童将购进的n份全部售出，利润为2()()nnbabnAnK。在两种情况下将利润与需求概率相乘并求和，就得到报童每天的平均利润()Vn：2210()[()()]()[()]()nrrnnnVnbcrAncnfrbnAnfrKK虽然报纸的需求量是离散的，但由于需求量很大，可以将它看作连续随机变量，所以报童每天的平均利润可以写作：220()[()()]()d[()]()dnnnnVnbcrAncnprrbnAnprrKK其中()pr是正态分布2(,)N的概率密度。为了求()Vn的最大，对()Vn求导数，得0022()[(1)]()[(1)]()d[(1)]()[(1)]()d22[(1)]()d[(1)]()dnnnnnnnnVnbAnpnAcprrbAnpnbAprrKKKKnnAcprrbAprrKK当()Vn取得最大值时，()0Vn，得到02(1)()d2()d(1)nnnbAprrKnprrAcK……(1)从题中可以看到200050，所以有0()d()dnnprrprr，又()d1()dnnprrprr，(1)式可化为2(1)()dnnbAKprrbc……(2)又因为002222()[(1)]()()d[(1)]()()d2()d20nnnnVnAcpnAprrbApnAprrKKKKcbAprrKcbAk所以当(2)式成立时得到的确是每天平均利润最大的最佳购进量n。利用MATLAB求解，输入m文件如下：mu=2000;sita=50;A=0.5;K=50000;b=0.5;c=0.35;n(:,1)=100;fork=1:50q=(b-A*(1-2*n(:,k)/K))/(b-c);n(:,k+1)=norminv(q,mu,sita);ifabs(n(:,k+1)-n(:,k))1e-1breakendendN=n(:,k+1)解得n=1968.2核算得n=1968时，报童的平均获利最大。结果：报童购进量为1968份时，其平均获利最大，即其盈利的期望最大。下面计算最大平均获利：2219680196819680019680.64004019681968(1968)[()(1968)1968]()d[1968(1968)]()d(0.15256.47)()d38.73()d(0,1)[0.15()256.47]()d38.73(VbcrAcprrbAprrKKrrprrprrxNxpxxpx令00.640.6400400.64)d[7.543.53]()d38.73()dxxpxxpxx将正态分布2(,)N的概率密度()pr换为(0,1)N的概率密度()x则有：220.64400.640.640.6440400.640.640.64240400.640.240(1968)[7.543.53]()d38.73()d7.5()d43.53()d38.73()d17.5d43.53()d38.73()d217.52xxVxxxxxxxxxxxxxexxxxxe640.64400.6443.53()d38.73()dxxxx利用MATLAB计算得出functionf=exe1007f(x)f=x*1/(2*pi)^0.5*exp(-x^2/2);functiong=exe1007g(x)g=1/(2*pi)^0.5*exp(-x^2/2);V1=-7.5*1/(2*pi)^0.5*(exp(-0.64^2/2)-exp(-40^2/2))+43.53*(normcdf(-0.64,0,1)-normcdf(-40,0,1))+38.73*(1-normcdf(-0.64,0,1))symsxV11=7.5*int(exe1007f(x),-40,-0.64)+43.53*int(exe1007g(x),-40,-0.64)+38.73*int(exe1007g(x),-0.64,inf)两种方法解得：V=V1=V11=37.5452。可以看到报童进购了1968份报纸，其最大平均获利才只有37.5452元。若报童刚好卖完，则实际获利最大为2()()nnbabnAnK，代入数据得38.73元。平均获利是指报童获利的期望，而不是报童的实际获利，报童的实际获利是随机的，其均值为V。结果分析（1）由于批发价为(1)naAK，所以n不能过大，否则报纸发行商将亏本；报纸发行商不仅需要限制n，还需要限制K的大小，K太小，会使得报童在购进较少n时，批发价就小于回收价格，造成亏本。下面考察K对结果的影响：利用MATLAB求解不同K时，报童的最优购进量如下：K最优购进量n50000196840000197820000202210000NaN可以看到，当K减小时，最优购进量增大，这是因为成倍减小，使得报童提高购进量所承担的亏本风险降低。当K减小到10000一下时，输出值为NaN，此时不存在最优值，n越大越好，报童不会亏本，但发行商亏本，所以K值不能过小。（2）n很接近，这很容易理解，当n时，报童报纸卖不出去的可能性较小，报童应尽量提高n使其能买更多的报纸；当n时，n越大，卖不出去的可能性就越大，所以n接近。下面考察对结果的影响：最优购进量n20001968100094250042050005022可以看到n很接近，与预测相符合，且越大，n就越大，这是因为此时n增大，购进成本(1)naAK减小，报童购进报纸承担的风险降低。可以看到在为5000时，n。（3）报纸需求量的方差是报童盈亏的风险因子，越大，需求量就越难预测，造成报童的风险。下面考察对结果的影响：最优购进量n5019682019871019931001935可以看到，当越大时，报童越本能多购进，因为需求量