MATLAB数学建模与仿真（第2版·微课视频版）程序代码 综合实例（十）

【例10-49】旱灾土地总面积问题。下表10-6是从1990年至2010年全国因干旱而受灾的土地总面积（单位：千公顷）数。（数据来源于全国统计局官网）表10-6各个年份全国受灾土地总面积年份受灾面积年份受灾面积年份受灾面积199018175199733516200417253199124914199814236200516028199232981199930156200620738199321097200040541200729386199430423200138472200812137199523455200222124200929259199620152200324852201013259解决以下问题：（1）计算所给样本的均值与标准差；（2）检验在显著水平为0.05的情况下，全国每年因干旱而受灾的土地总面积是否服从正态分布；（3）如果服从正态分布，用极大似然估计法对未知参数μ和σ作出估计；（4）若年受旱灾总面积大于35000千公顷即为重灾年，根据估计出的μ值和σ值，计算当年为重灾年的概率。分析问题：这是一个样本均值和标准差的计算以及正态性检验和计算的一系列问题。对于此类问题可以应用数学软件MATLAB进行处理，应用MATLAB可以很容易的计算出均值及标准差，此外，采用Jarque-Beran检验即可知道其是否服从正态分布，并估计出总体的均值μ和标准差σ。解决问题：下面计算样本的均值和标准差>>X=[181752491732981210973042323455201523351614236301564054138472221242485217253160282073829386121372925913259];>>[h,stats]=cdfplot(X)h=Line-属性:Color:[00.44700.7410]LineStyle:'-'LineWidth:0.5000Marker:'none'MarkerSize:6MarkerFaceColor:'none'XData:[1×44double]YData:[1×44double]ZData:[1×0double]显示所有属性stats=包含以下字段的struct:min:12137max:40541mean:2.4436e+04median:23455std:8.1234e+03图10-17样本经验分布函数图从输出结果可看出，样本的最小值为12137，最大值为40541，中值为23455，均值为24436，标准差为8123.4。验证下面检验其是否服从正态分布。MATLAB程序代码如下>>clearall;>>X=[181752491732981210973042323455201523351614236301564054138472221242485217253160282073829386121372925913259];>>normplot(X);>>[h,P,Jbstat,CV]=jbtest(X,0.05)>>title('正态概率图')>>xlabel('数据');>>ylable('概率')运行程序后，输出如下：h=0P=0.4574Jbstat=0.9230CV=3.8801图10-18正态概率图由输出结果h=0且Jbstat<CV可得出结论，在置信度α=0.05下，受灾面积（原始数据）服从正态分布，且在正态概率图中，各点均落在直线两侧，也可说明这一结论是成立的。再用极大似然估计法对未知参数μ和σ作出估计：MATLAB程序代码如下>>clearall;>>X=[181752491732981210973042323455201523351614236301564054138472221242485217253160282073829386121372925913259];>>phat=mle(X,'distribution','norm','alpha',0.05)运行程序后，输出如下phat=1.0e+04*2.44360.7928即受灾面积的μ估计值为24436，σ估计值为7928。最后根据估计出的μ值和σ值，计算出每年的受灾面积大于35000千公顷的概率：MATLAB程序代码如下>>clearall;p=normspec([35000inf],24436,7928)运行程序后，输出如下p=0.0913图10-19密度函数图根据输出结果可知，为重灾年的概率为0.0913.学习总结：通过对1990年至2010年全国因干旱而受灾的土地总面积的分析，我们得出这些数据服从正态分布。运用MATLAB程序，得出年均受灾土地总面积为24436公顷，根据图表可清晰地看出每年受灾总面积的分布状况，可以根据对这些数据的具体分布，采取相应的措施，从而最大程度的减小受灾。经过此次对实际问题的解决，让我们共同认识到概率统计的知识在我们的生活中无处不在，概率论在我们学习和生活中的应用也给人们带来了极大地便利。在对数据处理的过程中，对于很多数学工具的应用，如MATLAB等数学软件可让数据处理变得更加简单，从而引导我们更深层次的去探讨数学问题。在小组合作中让我们体会到小组分工合作的重要性，让我们受益匪浅。【例10-50】已知x=0:0.1:1，y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2]，求用polyfit和polyval函数拟合出的曲线图。%在命令行窗口输入如下命令>>x=0:0.1:1;>>y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];>>A=polyfit(x,y,2);>>z=polyval(A,x);>>plot(x,y,'r*',x,z,'b')图10-13曲线拟合的函数图【例10-51】以一简单数据组来说明什么是线性回归。假设有一组数据型态为y=y(x)，其中x={0,1,2,3,4,5},y={0,20,60,68,77,110}。如果要以一个最简单的方程式来近似这组数据，则用一阶的线性方程式最为适合。从polyfityy值间误差平方的总合（见图10-2。>>x=[012345];>>y=[020606877110];>>coef=polyfit(x,y,1);%coef代表线性回归的二个输出值>>a0=coef(1);a1=coef(2);>>ybest=a0*x+a1;%由线性回归产生的一阶方程式>>sum_sq=sum((y-ybest).^2);%误差平方总合为356.82>>axis([-1,6,-20,120])>>plot(x,ybest,x,y,'o'),title('Linearregressionestimate'),grid图10-14曲线拟合的函数图【例10-52】polyval计算在X中任意元素处的多项式p的估值。对多项式p(x)=1+2*x+3*x^2，计算在x=5、7、9的值。>>x=[5,7,9];>>p=[3,2,1];>>polyval(p,x)%结果为ans=86 162 262【例10-53】如何估计多元线性回归的系数，以MATLAB自带的数据为样本，示例代码如下。loadcarsmall%此数据样本MATLAB自带x1=Weight;%取这3个变量作为拟合对象，x1、x2自变量，y应变量x2=Horsepower;y=MPG;%用相互作用项计算线性模型的回归系数X=[ones(size(x1))x1x2x1.*x2];b=regress(y,X);%绘制数据和模型scatter3(x1,x2,y,'filled')holdonx1fit=min(x1):100:max(x1);x2fit=min(x2):10:max(x2);[X1FIT,X2FIT]=meshgrid(x1fit,x2fit);YFIT=b(1)+b(2)*X1FIT+b(3)*X2FIT+b(4)*X1FIT.*X2FIT;mesh(X1FIT,X2FIT,YFIT)xlabel('Weight')ylabel('Horsepower')zlabel('MPG')view(50,10)图10-15估计多元线性回归的系数图【例10-54】已知，x=[1.13891.06220.98220.9340.92510.9158]; y=[0.0315.0315.0519.97其中拟合函数为y=-k*ln(x+a)-b，求利用nlinfit函数进行非线性参数拟合。x=[1.13891.06220.98220.9340.92510.9158];y=[0.0315.0315.0519.9730.3];myfunc=inline('-beta(1)*log(x+beta(2))-beta(3)','beta','x');beta=nlinfit(x,y,myfunc,[000]);k=beta(1),a=beta(2),b=beta(3)%testthemodelxx=min(x):max(x);yy=-k*log(x+a)-b;plot(x,y,'o',x,yy,'r')图10-16nlinfit非线性参数拟合图【例10-55】已知，x=2:10;y=8*sin(x).*exp(x)-12./log(x)。现在已经有