MATLAB使用教程5應(yīng)用舉例課件

（五）MATLAB應(yīng)用舉例例5.12005高教社杯全國大學(xué)生數(shù)學(xué)建模競賽題目C雨量預(yù)報方法的評價

雨量預(yù)報對農(nóng)業(yè)生產(chǎn)和城市工作和生活有重要作用，但準確、及時地對雨量作出預(yù)報是一個十分困難的問題，廣受世界各國關(guān)注。我國某地氣象臺和氣象研究所正在研究6小時雨量預(yù)報方法，即每天晚上20點預(yù)報從21點開始的4個時段（21點至次日3點，次日3點至9點，9點至15點，15點至21點）在某些位置的雨量，這些位置位于東經(jīng)120度、北緯32度附近的53×47的等距網(wǎng)格點上。同時設(shè)立91個觀測站點實測這些時段的實際雨量，由于各種條件的限制，站點的設(shè)置是不均勻的。12/31/20220黃建華制作（五）MATLAB應(yīng)用舉例例5.12005高教社例5.1雨量預(yù)報方法的評價氣象部門希望建立一種科學(xué)評價預(yù)報方法好壞的數(shù)學(xué)模型與方法。氣象部門提供了41天的用兩種不同方法的預(yù)報數(shù)據(jù)和相應(yīng)的實測數(shù)據(jù)。預(yù)報數(shù)據(jù)在文件夾FORECAST中，實測數(shù)據(jù)在文件夾MEASURING中，其中的文件都可以用Windows系統(tǒng)的“寫字板”程序打開閱讀。FORECAST中的文件lon.dat和lat.dat分別包含網(wǎng)格點的經(jīng)緯度，其余文件名為<f日期i>_dis1和<f日期i>_dis2，例如f6181_dis1中包含2002年6月18日晚上20點采用第一種方法預(yù)報的第一時段數(shù)據(jù)（其2491個數(shù)據(jù)為該時段各網(wǎng)格點的雨量），而f6183_dis2中包含2002年6月18日晚上20點采用第二種方法預(yù)報的第三時段數(shù)據(jù)。12/31/20221黃建華制作例5.1雨量預(yù)報方法的評價氣象部門希望建立一種科學(xué)評價預(yù)報例5.1雨量預(yù)報方法的評價MEASURING中包含了41個名為<日期>.SIX的文件，如020618.SIX表示2002年6月18日晚上21點開始的連續(xù)4個時段各站點的實測數(shù)據(jù)（雨量），這些文件的數(shù)據(jù)格式是：站號緯度經(jīng)度第1段第2段第3段第4段5813832.9833118.51670.00000.200010.10003.10005813933.3000118.85000.00000.00004.60007.40005814133.6667119.26670.00000.00001.10001.40005814333.8000119.80000.00000.00000.00001.80005814633.4833119.81670.00000.00001.50001.9000……雨量用毫米做單位，小于0.1毫米視為無雨。12/31/20222黃建華制作例5.1雨量預(yù)報方法的評價MEASURING例5.1雨量預(yù)報方法的評價（1）請建立數(shù)學(xué)模型來評價兩種6小時雨量預(yù)報方法的準確性；（2）氣象部門將6小時降雨量分為6等：0.1—2.5毫米為小雨，2.6—6毫米為中雨，6.1—12毫米為大雨，12.1—25毫米為暴雨，25.1—60毫米為大暴雨，大于60.1毫米為特大暴雨。若按此分級向公眾預(yù)報，如何在評價方法中考慮公眾的感受？（注：本題數(shù)據(jù)位于壓縮文件C2005Data.rar中,可從/mcm05/problems2005c.asp下載）

12/31/20223黃建華制作例5.1雨量預(yù)報方法的評價（1）請建立數(shù)學(xué)模型來評價兩種6例5.1雨量預(yù)報方法的評價模型的分析：本題的關(guān)鍵主要是采用Matlab軟件對所提供的數(shù)據(jù)進行分析并得到主要結(jié)論。對于問題一，采用load命令和循環(huán)結(jié)構(gòu)實現(xiàn)數(shù)據(jù)文件的載入，再從實測數(shù)據(jù)文件中提出實測點位置數(shù)據(jù)，并且依次從預(yù)測數(shù)據(jù)文件中通過曲面擬合命令griddata得到相應(yīng)日期、時段、方法下的對應(yīng)位置上的預(yù)測估計值。然后分別計算兩種預(yù)報方法下的實測值和預(yù)測值的偏差并且求對應(yīng)的總偏差平方和，根據(jù)兩個總偏差平方和的大小來得出兩種方法的優(yōu)劣比較，通過Matlab程序的實現(xiàn)得到結(jié)果為：第一種方法比第二種方法好。12/31/20224黃建華制作例5.1雨量預(yù)報方法的評價模型的分析：12/26/2022例5.1雨量預(yù)報方法的評價對于問題二，我們認為公眾的感受主要體現(xiàn)在預(yù)測與實際之間偏差的大小程度，然而，實際測量值未知，因此我們考慮在分級預(yù)報中加入準確概率的方式來實現(xiàn)公眾滿意度的提高。我們認為雨量的實測值應(yīng)該在預(yù)測值點處服從正態(tài)分布，通過合理的假設(shè)和推導(dǎo)，我們得到：正態(tài)分布的均值可以取為雨量預(yù)測值，不同雨量分級區(qū)間上的方差可以近似取為對應(yīng)區(qū)間上的總偏差平方和的平均值。運用Matlab軟件編程，可以實現(xiàn)對每一個預(yù)報數(shù)據(jù)的預(yù)報內(nèi)容的改變，并且通過幾個不同的數(shù)據(jù)體現(xiàn)程序運行所得到的結(jié)果。最后，對模型的缺點進行了討論和改進。

12/31/20225黃建華制作例5.1雨量預(yù)報方法的評價對于問題二，我們認為公眾的感受主例5.1雨量預(yù)報方法的評價模型的假設(shè)假設(shè)測量雨量的工具正常不受任何因素的影響，且所得數(shù)據(jù)真實可靠；測量點所在位置的等距網(wǎng)格點均視為質(zhì)點；符號說明x(i,j)表示第j種方法下的第i個實測點處的雨量實測值y(i,j)表示第j種方法下的第i個實測點處的雨量預(yù)測估計值i=1,2,3…….,91*164,j=1,2e(i,j)表示在第j種方法下的第i個實測點處的偏差Y表示實際值的隨機變量y*表示預(yù)報值ST2(j)表示在第j種方法下的所有實測點處的總偏差平方和12/31/20226黃建華制作例5.1雨量預(yù)報方法的評價模型的假設(shè)12/26/20226例5.1雨量預(yù)報方法的評價問題一模型分析：考察下載的lon.dat、lat.dat、020<日期>.SIX及f<日期><i>_dis<j>等數(shù)據(jù)文件。我們可以看到雨量預(yù)報的網(wǎng)格點有53*47個。而雨量實測點只有91個，且分布不均勻，因此我們不能直接引用數(shù)據(jù)文件進行實測數(shù)據(jù)與預(yù)測數(shù)據(jù)之間偏差的計算。顯然首先要對所提供的數(shù)據(jù)文件進行處理，由于所引用的文件比較多，而且具有一定的規(guī)律性，所以數(shù)據(jù)文件的載入可用MATLAB命令load和循環(huán)結(jié)構(gòu)來實現(xiàn)。當(dāng)有了具體的實測點位置對應(yīng)的實測數(shù)據(jù)yi及預(yù)測數(shù)據(jù)yi*之后,我們可以求出每種方法下所有實測點位置上的總偏差平方和。然后比較兩個總偏差平方和的大小，就可以判斷兩種方法的優(yōu)劣性。

12/31/20227黃建華制作例5.1雨量預(yù)報方法的評價問題一模型分析：12/26/20例5.1雨量預(yù)報方法的評價問題一模型建立：

通過Matlab程序編程計算，具體過程如下：（1）使用MATLAB命令load和循環(huán)結(jié)構(gòu)將lot.dat、lat.dat、f6181_dis1…、f7304_dis2、020618.six…、020730.six等文件輸入MATLAB程序中備用，其中020618.six…、020730.six等文件系統(tǒng)自動在前面添上X。（2）從X020618文件中取第二、第三列作為實測點位置。其中第一列對應(yīng)實測點的緯度，第二列對應(yīng)實測點的經(jīng)度。分別用lat、lon中的數(shù)據(jù)作為網(wǎng)格點的橫坐標和縱坐標。依次取f6181_dis1……,f7304_dis1,f6181_dis2,……,f7304_dis2文件對應(yīng)點數(shù)據(jù)作為豎坐標，張成一個預(yù)報值數(shù)據(jù)曲面，用曲面擬合命令griddata得出對應(yīng)實測點處的預(yù)報估計值，得到一個12/31/20228黃建華制作例5.1雨量預(yù)報方法的評價問題一模型建立：12/26/2例5.1雨量預(yù)報方法的評價91*328維的對應(yīng)實測點處的預(yù)報估計值矩陣yczjz。（3）依次從X020618，…….X020730文件中取第四到第七列組成所有日期，所有時段下的實測數(shù)據(jù)矩陣sczjz。（4）將sczjz分別與yczjz前164列、后164列對應(yīng)元素相減，得出第一種方法下和第二種方法下所有實測點處的偏差，組成矩陣f1pcz和f2pcz。（5）分別求矩陣f1pcz和f2pcz所有元素的平方之和，賦予變量pcpfh(1)和pcpfh(2)。（6）比較變量pcpfh(1)、pcpfh(2)的大小，得出兩種方法優(yōu)劣的比較。各過程的MATLAB程序?qū)崿F(xiàn)，可分別由M-文件sjzr.m、qsj.m、pcbj.m運行得到。

12/31/20229黃建華制作例5.1雨量預(yù)報方法的評價91*328維的對應(yīng)實測點處的預(yù)例5.1雨量預(yù)報方法的評價程序sjzr.m：

rqjz=[618:628,701:730];%產(chǎn)生1*41維日期行向量loadlat.DAT;%載入53*47維緯度矩陣loadlon.DAT;%載入53*47維經(jīng)度矩陣forrq=rqjz%用日期行向量作為循環(huán)向量load(['020',int2str(rq),'.SIX'])%載入53*47維實測值矩陣,系統(tǒng)自動在文件名前加大寫字母Xforsd=1:4%指定時段循環(huán)向量forff=1:2%指定方法循環(huán)向量load(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)])%載入53*47維預(yù)測值矩陣，共328個endendend

12/31/202210黃建華制作例5.1雨量預(yù)報方法的評價程序sjzr.m：12/26/例5.1雨量預(yù)報方法的評價將文件全部拷貝在默認的work文件加下，運行程序sjzr.m：

12/31/202211黃建華制作例5.1雨量預(yù)報方法的評價將文件全部拷貝在默認的work文例5.1雨量預(yù)報方法的評價12/31/202212黃建華制作例5.1雨量預(yù)報方法的評價12/26/202212黃建華制例5.1雨量預(yù)報方法的評價程序huatu.m%畫一個時點預(yù)測曲面圖與實測散點圖比較clf;%預(yù)先清除別的圖象Z=eval(['f',int2str(618),int2str(1),'_dis',int2str(1)]);%取出某時段預(yù)報值作為豎坐標mesh(lat,lon,Z)%張成預(yù)測曲面axis([27,36,117,125,0,0.4])shadingflatxlabel('緯度');ylabel('經(jīng)度');title('圖象比較');holdonX1=X020618(:,2);Y1=X020618(:,3);Z1=X020618(:,4);%取出實測點數(shù)據(jù)plot3(X1,Y1,Z1,'r*')%畫出實測點數(shù)據(jù)散點圖12/31/202213黃建華制作例5.1雨量預(yù)報方法的評價程序huatu.m%畫一例5.1雨量預(yù)報方法的評價12/31/202214黃建華制作例5.1雨量預(yù)報方法的評價12/26/202214黃建華制例5.1雨量預(yù)報方法的評價程序qsj.mycr=1;%指定預(yù)測值矩陣列標yczjz=zeros(91,328);%預(yù)先指定預(yù)測值矩陣為全零矩陣forff=1:2%指定方法循環(huán)變量forrq=rqjz%指定日期循環(huán)變量forsd=1:4%指定時段循環(huán)變量z0=eval(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)]);%依次取出雨量預(yù)測值數(shù)據(jù)yczjz(:,ycr)=griddata(lat,lon,z0,X020618(:,2),X020618(:,3));%在張成的預(yù)測數(shù)據(jù)曲面上擬合對應(yīng)實測點的雨量預(yù)測值，依次放入預(yù)測值矩陣對應(yīng)列ycr=ycr+1;%預(yù)測值矩陣列標向后一列end12/31/202215黃建華制作例5.1雨量預(yù)報方法的評價程序qsj.m12/26/202例5.1雨量預(yù)報方法的評價endendscr=1;%指定實測值矩陣列標sczjz=zeros(91,164);%預(yù)先指定實測值矩陣為全零矩陣forrq=rqjz%指定日期循環(huán)變量sczjz(:,scr:scr+3)=eval(['X020',int2str(rq),'(:,4:7)']);%依次取出雨量實測值數(shù)據(jù)放入實測值矩陣scr=scr+4;%實測值矩陣列標向后4列end

12/31/202216黃建華制作例5.1雨量預(yù)報方法的評價end12/26/2例5.1雨量預(yù)報方法的評價程序pcbj.mff1=yczjz(:,1:164);%取出預(yù)測值矩陣的前164列ff2=yczjz(:,165:328);%取出預(yù)測值矩陣的后164列f1pcz=sczjz-ff1;%計算方法1下的偏差f2pcz=sczjz-ff2;%計算方法2下的偏差pcpfh=[0,0];%預(yù)先產(chǎn)生1*2維偏差平方和零矩陣forpch=1:91%指定偏差行循環(huán)變量forpcr=1:164%指定偏差列循環(huán)變量pcpfh(1)=pcpfh(1)+f1pcz(pch,pcr).^2;%計算方法1偏差平方和pcpfh(2)=pcpfh(2)+f2pcz(pch,pcr).^2;%計算方法2偏差平方和12/31/202217黃建華制作例5.1雨量預(yù)報方法的評價程序pcbj.m12/26/20例5.1雨量預(yù)報方法的評價endendpcpfhifpcpfh(1)==pcpfh(2)%比較兩種方法下的偏差平方和display('一樣好')elseifpcpfh(1)<pcpfh(2)display('第一種方法好')elsedisplay('第二種方法好')end首先將數(shù)據(jù)文件減壓縮在matlab程序的work文件夾下，然后依次運行M-文件sjzr.m、qsj.m、pcbj.m，得到如下結(jié)論：偏差平方和：226730243880第一種方法比第二種方法好。12/31/202218黃建華制作例5.1雨量預(yù)報方法的評價end12/26/202例5.1雨量預(yù)報方法的評價問題二模型分析：在對問題二的分析中，我們只對方法一相應(yīng)的數(shù)據(jù)進行分析。在現(xiàn)實生活中，公眾對天氣預(yù)報的感受主要來自于對雨量預(yù)報等級和雨量實際等級之間的偏差所帶來的不滿，而不是對雨量等級本身的不滿，因此我們認為，考慮公眾的感受就是要盡量減少所預(yù)報等級和實際等級之間偏差的程度。但在預(yù)報過程中，我們不可能引用來自于還未知的實際數(shù)據(jù)，我們只能從評價方法中提供的預(yù)報數(shù)據(jù)及以前的歷史數(shù)據(jù)來提高預(yù)報的準確性及公眾的感受度。為此我們考慮在分級預(yù)報中引入概率。比如：預(yù)報“中雨”時，我們可以采取預(yù)報方式為“中雨，概率為：55%”等。12/31/202219黃建華制作例5.1雨量預(yù)報方法的評價問題二模型分析：12/26/20例5.1雨量預(yù)報方法的評價一般的，雨量實際值是在預(yù)報值附近振動的，而且越接近預(yù)報值的實際值的概率越高，離預(yù)報值越遠的實際值的概率越低，因此我們可以認為：實際值服從以預(yù)報值為中心的正態(tài)分布。假設(shè)實際值隨機變量為Y，預(yù)報值y*則Y~N(y*,σ2)，其中σ2待定?，F(xiàn)假設(shè)對實際值Y進行n次檢測得到樣本為Yi，I=1,2,…..,n,對應(yīng)樣本值為yi，顯然，由概率統(tǒng)計知識可知Yi~N(y*,σ2),且是y*的無偏估計，是σ2的無偏估計。即可用代替y*，s2代替σ2。12/31/202220黃建華制作例5.1雨量預(yù)報方法的評價一般的，雨量實際值是在預(yù)報例5.1雨量預(yù)報方法的評價則：其中ei為每次預(yù)報的偏差。在此，我們顯然不能對同一位置進行n次偏差檢測。但是可以用在同一分級下的實測值與預(yù)測值之間偏差近似代替ei并求σ2值。然后再用各分級下的σ2值求給定預(yù)報值所在對應(yīng)分級下的預(yù)報準確概率P。

12/31/202221黃建華制作例5.1雨量預(yù)報方法的評價則：12/26/202221例5.1雨量預(yù)報方法的評價問題二模型建立：

由上述分析可知，若將無雨、小雨、中雨、大雨、暴雨、大暴雨、特大暴雨等7個等級所對應(yīng)的降雨量區(qū)間記為Ij=[aj,bj],則當(dāng)y*∈Ij時，在第j分級下的均方差為σ2j，Y~N(y*,σ2j),而且Y的實際值y∈Ij的概率為：P=預(yù)報方式為：等級，概率為：P12/31/202222黃建華制作例5.1雨量預(yù)報方法的評價問題二模型建立：12/26/2例5.1雨量預(yù)報方法的評價程序pcbj.mff1=yczjz(:,1:164);%取出預(yù)測值矩陣的前164列ff2=yczjz(:,165:328);%取出預(yù)測值矩陣的后164列f1pcz=sczjz-ff1;%計算方法1下的偏差f2pcz=sczjz-ff2;%計算方法2下的偏差pcpfh=[0,0];%預(yù)先產(chǎn)生1*2維偏差平方和零矩陣forpch=1:91%指定偏差行循環(huán)變量forpcr=1:164%指定偏差列循環(huán)變量pcpfh(1)=pcpfh(1)+f1pcz(pch,pcr).^2;%計算方法1偏差平方和pcpfh(2)=pcpfh(2)+f2pcz(pch,pcr).^2;%計算方法2偏差平方和12/31/202223黃建華制作例5.1雨量預(yù)報方法的評價程序pcbj.m12/26/20例5.1雨量預(yù)報方法的評價程序執(zhí)行過程描述如下：（1）用if-else-end語句和for循環(huán)對問題一中所得的偏差矩陣進行分級求偏差平方和，并累計個數(shù)。（2）計算各分級下的偏差平方和與長度減一的比值，分別作為對應(yīng)等級下的正態(tài)分布的均方差的近似值。（3）輸入預(yù)測值，利用求正態(tài)分布分布函數(shù)命令normcdf及上步中σ值計算對應(yīng)分級中的預(yù)報準確概率。相應(yīng)程序見M-文件flpcjz.m、yb.m.12/31/202224黃建華制作例5.1雨量預(yù)報方法的評價程序執(zhí)行過程描述如下：12/26例5.1雨量預(yù)報方法的評價程序flpcjz.mwypcf=0;%給‘無雨’偏差平方和賦初值xypcf=0;%給‘小雨’偏差平方和賦初值zypcf=0;%給‘中雨’偏差平方和賦初值dypcf=0;%給‘大雨’偏差平方和賦初值bypcf=0;%給‘暴雨’偏差平方和賦初值dbypcf=0;%給‘大暴雨’偏差平方和賦初值tdypcf=0;%給‘特大暴雨’偏差平方和賦初值wycd=0;%給‘無雨’偏差個數(shù)賦初值xycd=0;%給‘小雨’偏差個數(shù)賦初值zycd=0;%給‘中雨’偏差個數(shù)賦初值dycd=0;%給‘大雨’偏差個數(shù)賦初值bycd=0;%給‘暴雨’偏差個數(shù)賦初值12/31/202225黃建華制作例5.1雨量預(yù)報方法的評價程序flpcjz.m12/26/例5.1雨量預(yù)報方法的評價dbycd=0;%給‘大暴雨’偏差個數(shù)賦初值tdycd=0;%給‘特大暴雨’偏差個數(shù)賦初值forcol=1:91%用預(yù)測值矩陣行數(shù)作為循環(huán)變量forrow=1:164,%用預(yù)測值矩陣列數(shù)作為循環(huán)變量ifyczjz(col,row)<=0.1%判斷預(yù)測值矩陣對應(yīng)位置處值的分類wypcf=wypcf+f1pcz(col,row).^2;wycd=wycd+1;%按分類累加偏差，并且對應(yīng)個數(shù)加1elseifyczjz(col,row)>0.1&yczjz(col,row)<=2.5,xypcf=xypcf+f1pcz(col,row).^2;xycd=xycd+1;elseifyczjz(col,row)>2.5&yczjz(col,row)<=6,zypcf=zypcf+f1pcz(col,row).^2;zycd=zycd+1;12/31/202226黃建華制作例5.1雨量預(yù)報方法的評價dbycd=0;例5.1雨量預(yù)報方法的評價elseifyczjz(col,row)>6&yczjz(col,row)<=12,dypcf=dypcf+f1pcz(col,row).^2;dycd=dycd+1;elseifyczjz(col,row)>12&yczjz(col,row)<25,bypcf=bypcf+f1pcz(col,row).^2;bycd=bycd+1;elseifyczjz(col,row)>25&yczjz(col,row)<=60,dbypcf=dbypcf+f1pcz(col,row).^2;dbycd=dbycd+1;elsetdypcf=tdypcf+f1pcz(col,row).^2;tdycd=tdycd+1;endendend12/31/202227黃建華制作例5.1雨量預(yù)報方法的評價elseifyczjz(co例5.1雨量預(yù)報方法的評價程序yb.mwz1='無雨';wz2='小雨';wz3='中雨';wz4='大雨';wz5='暴雨';wz6='大暴雨';wz7='特大暴雨';ifycz<=0.1gl=normcdf(0.1,ycz,sgm(1));wz=wz1;%根據(jù)預(yù)報值計算分級和相應(yīng)概率elseifycz>0.1&ycz<=2.5wz=wz2;12/31/202228黃建華制作例5.1雨量預(yù)報方法的評價程序yb.m12/26/2022例5.1雨量預(yù)報方法的評價gl=normcdf(2.5,ycz,sgm(2))-normcdf(0.1,ycz,sgm(2));elseifycz>2.5&ycz<=6gl=normcdf(6,ycz,sgm(3))-normcdf(2.5,ycz,sgm(3));wz=wz3;elseifycz>6&ycz<=12gl=normcdf(12,ycz,sgm(4))-normcdf(6,ycz,sgm(4));wz=wz4;elseifycz>12&ycz<=25gl=normcdf(25,ycz,sgm(5))-normcdf(12,ycz,sgm(5));wz=wz5;elseifycz>25&ycz<=60gl=normcdf(60,ycz,sgm(6))-normcdf(25,ycz,sgm(6));wz=wz6;12/31/202229黃建華制作例5.1雨量預(yù)報方法的評價gl=normcd例5.1雨量預(yù)報方法的評價elsegl=1-normcdf(60,ycz,sgm(7));wz=wz7;enddisplay([eval('wz'),',概率為:']),display(gl)%輸出等級和概率

首先運行程序flpcjz.m，得到：sgm=[0.08070.60641.95762.74535.191513.5697105.4409]然后：輸入預(yù)報值：ycz=2.2結(jié)果為：小雨，概率為：0.6893ycz=0.2結(jié)果為：小雨，概率為：0.5654ycz=4.5結(jié)果為：中雨，概率為：0.6248ycz=15.6結(jié)果為：暴雨，概率為：0.720912/31/202230黃建華制作例5.1雨量預(yù)報方法的評價else12/26/2022例5.1雨量預(yù)報方法的評價模型的進一步討論：

(1)在對matlab程序進行調(diào)試時，我們觀察發(fā)現(xiàn)矩陣中出現(xiàn)負值，但顯然對天氣預(yù)報數(shù)據(jù)來說，降雨量是不可能存在負值的，這個缺陷是由插值方法的不足所引起的，在程序中我們對此做了一定的改進。見程序qsjxg.m，通過設(shè)置一個條件結(jié)構(gòu)，當(dāng)預(yù)測雨量插值為正時，保留原值不變；當(dāng)預(yù)測雨量插值為負時，我們?nèi)☆A(yù)測雨量插值為0。經(jīng)過此改進，我們得到相應(yīng)的問題一的答案為：偏差平方和：226690243840第一種方法比第二種方法好。與原答案（偏差平方和：226730243880第一種方法比第二種方法好）相比較，可見插值誤差對結(jié)果影響不大。12/31/202231黃建華制作例5.1雨量預(yù)報方法的評價模型的進一步討論：12/26/例5.1雨量預(yù)報方法的評價(2)觀察問題二的演示答案，我們發(fā)現(xiàn)有些預(yù)報的概率有接近于50%的情況存在，表明有另外一種雨量等級概率比較大的可能。為了增加預(yù)報等級的可信度，不妨可以采取如下預(yù)報形式：大到暴雨，大雨概率：55%，暴雨概率40%。只是計算量有所增加，程序復(fù)雜程度有所增加。12/31/202232黃建華制作例5.1雨量預(yù)報方法的評價(2)觀察問題二的演示答案例5.1雨量預(yù)報方法的評價程序qsjxg.mycr=1;%指定預(yù)測值矩陣列標yczjz=zeros(91,328);%指定預(yù)測值矩陣為全零矩陣forff=1:2%指定方法循環(huán)變量forrq=rqjz%指定日期循環(huán)變量forsd=1:4%指定時段循環(huán)變量z0=eval(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)]);%依次取出雨量預(yù)測值數(shù)據(jù)yczjz(:,ycr)=griddata(lat,lon,z0,X020618(:,2),X020618(:,3));%在張成的預(yù)測數(shù)據(jù)曲面上擬合對應(yīng)實測點的雨量預(yù)測值，依次放入預(yù)測值矩陣對應(yīng)列forhb=1:91%指定行標循環(huán)變量ifyczjz(hb,ycr)>=0%判斷預(yù)測擬合值符合條件12/31/202233黃建華制作例5.1雨量預(yù)報方法的評價程序qsjxg.m12/26/2例5.1雨量預(yù)報方法的評價yczjz(hb,ycr)=yczjz(hb,ycr);%符合條件，則值不變elseyczjz(bjl,ycr)=0;%不符合條件，則值為0endendycr=ycr+1;%預(yù)測值矩陣列標向后一列endendendscr=1;%指定實測值矩陣列標forrq=rqjz%指定日期循環(huán)變量sczjz(:,scr:scr+3)=eval(['X020',int2str(rq),'(:,4:7)']);%依次取出雨量實測值數(shù)據(jù)scr=scr+4;%實測值矩陣列標加4end12/31/202234黃建華制作例5.1雨量預(yù)報方法的評價yczjz(hb,ycr)=yc（五）MATLAB應(yīng)用舉例例5.22004全國大學(xué)生數(shù)學(xué)建模競賽題目飲酒駕車的數(shù)學(xué)模型飲酒駕車問題主要是分析駕駛員在喝過一定量的酒后，酒精在體內(nèi)被吸收后，血液中酒精含量上升，影響司機駕車，所以司機飲酒后需經(jīng)過一段時間后才能安全駕車，國家標準新規(guī)定，車輛駕駛?cè)藛T血液中的酒精含量大于或等于20毫克/百毫升，小于80毫克/百毫升為飲酒駕車，血液中酒精含量大于或等于80毫克/百毫升為醉酒駕車，司機大李在中午12點喝下一瓶啤酒，6小時后檢查符合新標準，晚飯地其又喝了一瓶啤酒，他到凌晨2點駕車，被檢查時定為飲酒駕車，為什么喝相同量的酒，兩次結(jié)果不一樣？ 12/31/202235黃建華制作（五）MATLAB應(yīng)用舉例例5.22004全國大學(xué)例5.2飲酒駕車的數(shù)學(xué)模型討論問題：1、對大李碰到的情況做出合理解釋；2、在喝三瓶啤酒或半斤白酒后多長時間內(nèi)駕車會違反標準，喝酒時間長短不同情況會怎樣？3、分析當(dāng)司機喝酒后何時血液中的酒精含量最高；4、如果該司機想天天喝酒還能否開車；12/31/202236黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型討論問題：12/26/20223例5.2飲酒駕車的數(shù)學(xué)模型摘要本文解決的是一個司機安全駕車與飲酒的問題，目的是通過建立一個數(shù)學(xué)模型（結(jié)合新的國家駕駛員飲酒標準）分析司機如何適量飲酒不會影響正常的安全駕駛。根據(jù)一定合理的假設(shè)，建立人體內(nèi)酒精濃度隨時間變化的微分方程模型，并通過擬合曲線對數(shù)據(jù)進行分析。在不同飲酒方式下進行分類討論，得出體內(nèi)酒精濃度隨時間的變化函數(shù)。在討論過程中，我們得到兩個結(jié)論：在短時間喝酒形式下，達到最大值的時間為1.23小時，與喝酒量無關(guān)；在長時間喝酒形式下，喝酒結(jié)束時酒精含量最高。12/31/202237黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型摘要12/26/202237黃建華例5.2飲酒駕車的數(shù)學(xué)模型模型假設(shè)1、酒精從胃轉(zhuǎn)移到體液的速率與胃中的酒精濃度成正比。2、酒精從體液轉(zhuǎn)移到體外的速率與體液中的酒精濃度成正比。3、酒精從胃轉(zhuǎn)移到體液的過程中沒有損失。4、測量設(shè)備完善，不考慮不同因素所造成的誤差。5、酒精在體液中均勻分布。12/31/202238黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型模型假設(shè)12/26/202238黃例5.2飲酒駕車的數(shù)學(xué)模型符號說明k:酒精從體外進入胃的速率；f1(t):酒精從胃轉(zhuǎn)移到體液的速率；f2(t):酒精從體液轉(zhuǎn)移到體外的速率；X(t):胃里的酒精含量；Y(t)：體液中酒精含量；V0:體液的容積；K1:酒精從胃轉(zhuǎn)移到體液的轉(zhuǎn)移速率系數(shù)；K2:酒精從體液轉(zhuǎn)移到體外的轉(zhuǎn)移速率系數(shù)；C(t):體液中的酒精濃度。12/31/202239黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型符號說明12/26/202239黃例5.2飲酒駕車的數(shù)學(xué)模型D0：短時間喝酒情況下進入胃中的初始酒精量。T：較長時間喝酒所用的時間或達到濃度最大值所需時間。模型分析：假設(shè)酒精先以速率k0進入胃中，然后以速率f1(t)從胃進入體液，再以速率f2(t)從體液中排到體外。根據(jù)假設(shè)可以建立如圖所示的帶有吸收室的單房室系統(tǒng)，其中胃為吸收室，體液為中心室。12/31/202240黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型D0：短時間喝酒情況下進入胃中的初例5.2飲酒駕車的數(shù)學(xué)模型根據(jù)酒精從胃進入體液的速度f1(t)與胃中的酒精量成正比，速率系數(shù)為K1；酒精從血液中排出的速率f2(t)與血液中的酒精量y(t)成正比，速率系數(shù)為K2，可以建立方程如下：12/31/202241黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型根據(jù)酒精從胃進入體液的速度f1(t例5.2飲酒駕車的數(shù)學(xué)模型12/31/202242黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型12/26/202242黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型利用MATLAB的dsolve()函數(shù)，編輯：clear;symsxyy0k0k1k2;[x,y]=dsolve('Dx=k0-k1*x,Dy=k1*x-k2*y',…'x(0)=x0,y(0)=y0','t')pretty(simple(x)),pretty(simple(y))得到：12/31/202243黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型利用MATLAB的dsolve()例5.2飲酒駕車的數(shù)學(xué)模型12/31/202244黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型12/26/202244黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型模型討論：

當(dāng)酒是在較短時間內(nèi)喝時此時有x(0)=x0=D0，k0=0,y(0)=y0=0:又酒精濃度為酒精量與體液容積之比，即：12/31/202245黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型模型討論：12/26/202245例5.2飲酒駕車的數(shù)學(xué)模型顯然K1>>K2。當(dāng)t比較大時,可認為：利用數(shù)表：時間(小時)0.250.50.7511.522.533.544.55酒精含量153437.5414138.534342925.52520.5時間(小時)678910111213141516

酒精含量1917.51412.597.5653.53.52

12/31/202246黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型顯然K1>>K2。當(dāng)t比較大時,例5.2飲酒駕車的數(shù)學(xué)模型通過Matlab進行曲線擬合，編程：>>t=[1.522.533.544.55678910111213141516];y1=[82776868585150413835282518151210774];y2=log(y1);polyfit(t,y2,1);ans=-0.19404.7753可得：A=118.5459,k2=0.1940根據(jù)查閱資料可知：一瓶啤酒的酒精量一般為640ml，密度為810mg/ml酒精濃度為84.5%所以兩瓶啤酒的酒精總量由于體重為70kg,體重的65%左右，體液密度為1.05mg/ml，所以可得體液的總體積為12/31/202247黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型通過Matlab進行曲線擬合，編程例5.2飲酒駕車的數(shù)學(xué)模型可得短時間內(nèi)喝下兩瓶啤酒時血液中的酒精含量與時間的關(guān)系式如下:12/31/202248黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型可得短時間內(nèi)喝下兩瓶啤酒時血液中的例5.2飲酒駕車的數(shù)學(xué)模型用Matlab軟件畫出圖形，編程為：>>t=[0.250.50.7511.522.533.544.55678910…111213141516];y1=[30687582827768685851504138352825;a=118.5459;k1=2.114;k2=0.1940;x=0:20;y=a*(exp(-k2*t)-exp(-k1*t));plot(t,y1,'+',t,y,'r')12/31/202249黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型用Matlab軟件畫出圖形，編程為例5.2飲酒駕車的數(shù)學(xué)模型圖形為：12/31/202250黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型圖形為：12/26/202250黃例5.2飲酒駕車的數(shù)學(xué)模型問題一解答：假設(shè)大李第一次喝酒是在短時間內(nèi)喝的，根據(jù)所建立模型，所以可求得：

當(dāng)t=6時，可以求得c(t)=18.2778毫克/百毫升，小于國家規(guī)定的新標準，所以第一次遭遇檢查時沒有被認定為是飲酒駕駛，見圖：12/31/202251黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型問題一解答：12/26/20225例5.2飲酒駕車的數(shù)學(xué)模型圖形為：12/31/202252黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型圖形為：12/26/202252黃例5.2飲酒駕車的數(shù)學(xué)模型大李在吃晚飯時又喝了一瓶啤酒，喝酒后體內(nèi)酒精疊加.當(dāng)t=14時，c(t)=20.3618毫克/百毫升大于國家新規(guī)定的20mg/100ml,所以酒精超標：12/31/202253黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型大李在吃晚飯時又喝了一瓶啤酒，喝酒例5.2飲酒駕車的數(shù)學(xué)模型問題三解答：短時間內(nèi)喝酒時，根據(jù)所建立模型可知：當(dāng)c(t)的導(dǎo)數(shù)等于0時，可解得：因為T只與k1,k2有關(guān)，從表達式可知當(dāng)在較短時間內(nèi)喝酒時濃度達到最大值的時間與喝酒量無關(guān)。12/31/202254黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型問題三解答：12/26/20225例5.2飲酒駕車的數(shù)學(xué)模型其余解答請參考其他資料12/31/202255黃建華制作例5.2飲酒駕車的數(shù)學(xué)模型其余解答請參考其他資料12/26/（五）MATLAB應(yīng)用舉例例5.12005高教社杯全國大學(xué)生數(shù)學(xué)建模競賽題目C雨量預(yù)報方法的評價

雨量預(yù)報對農(nóng)業(yè)生產(chǎn)和城市工作和生活有重要作用，但準確、及時地對雨量作出預(yù)報是一個十分困難的問題，廣受世界各國關(guān)注。我國某地氣象臺和氣象研究所正在研究6小時雨量預(yù)報方法，即每天晚上20點預(yù)報從21點開始的4個時段（21點至次日3點，次日3點至9點，9點至15點，15點至21點）在某些位置的雨量，這些位置位于東經(jīng)120度、北緯32度附近的53×47的等距網(wǎng)格點上。同時設(shè)立91個觀測站點實測這些時段的實際雨量，由于各種條件的限制，站點的設(shè)置是不均勻的。12/31/202256黃建華制作（五）MATLAB應(yīng)用舉例例5.12005高教社例5.1雨量預(yù)報方法的評價氣象部門希望建立一種科學(xué)評價預(yù)報方法好壞的數(shù)學(xué)模型與方法。氣象部門提供了41天的用兩種不同方法的預(yù)報數(shù)據(jù)和相應(yīng)的實測數(shù)據(jù)。預(yù)報數(shù)據(jù)在文件夾FORECAST中，實測數(shù)據(jù)在文件夾MEASURING中，其中的文件都可以用Windows系統(tǒng)的“寫字板”程序打開閱讀。FORECAST中的文件lon.dat和lat.dat分別包含網(wǎng)格點的經(jīng)緯度，其余文件名為<f日期i>_dis1和<f日期i>_dis2，例如f6181_dis1中包含2002年6月18日晚上20點采用第一種方法預(yù)報的第一時段數(shù)據(jù)（其2491個數(shù)據(jù)為該時段各網(wǎng)格點的雨量），而f6183_dis2中包含2002年6月18日晚上20點采用第二種方法預(yù)報的第三時段數(shù)據(jù)。12/31/202257黃建華制作例5.1雨量預(yù)報方法的評價氣象部門希望建立一種科學(xué)評價預(yù)報例5.1雨量預(yù)報方法的評價MEASURING中包含了41個名為<日期>.SIX的文件，如020618.SIX表示2002年6月18日晚上21點開始的連續(xù)4個時段各站點的實測數(shù)據(jù)（雨量），這些文件的數(shù)據(jù)格式是：站號緯度經(jīng)度第1段第2段第3段第4段5813832.9833118.51670.00000.200010.10003.10005813933.3000118.85000.00000.00004.60007.40005814133.6667119.26670.00000.00001.10001.40005814333.8000119.80000.00000.00000.00001.80005814633.4833119.81670.00000.00001.50001.9000……雨量用毫米做單位，小于0.1毫米視為無雨。12/31/202258黃建華制作例5.1雨量預(yù)報方法的評價MEASURING例5.1雨量預(yù)報方法的評價（1）請建立數(shù)學(xué)模型來評價兩種6小時雨量預(yù)報方法的準確性；（2）氣象部門將6小時降雨量分為6等：0.1—2.5毫米為小雨，2.6—6毫米為中雨，6.1—12毫米為大雨，12.1—25毫米為暴雨，25.1—60毫米為大暴雨，大于60.1毫米為特大暴雨。若按此分級向公眾預(yù)報，如何在評價方法中考慮公眾的感受？（注：本題數(shù)據(jù)位于壓縮文件C2005Data.rar中,可從/mcm05/problems2005c.asp下載）

12/31/202259黃建華制作例5.1雨量預(yù)報方法的評價（1）請建立數(shù)學(xué)模型來評價兩種6例5.1雨量預(yù)報方法的評價模型的分析：本題的關(guān)鍵主要是采用Matlab軟件對所提供的數(shù)據(jù)進行分析并得到主要結(jié)論。對于問題一，采用load命令和循環(huán)結(jié)構(gòu)實現(xiàn)數(shù)據(jù)文件的載入，再從實測數(shù)據(jù)文件中提出實測點位置數(shù)據(jù)，并且依次從預(yù)測數(shù)據(jù)文件中通過曲面擬合命令griddata得到相應(yīng)日期、時段、方法下的對應(yīng)位置上的預(yù)測估計值。然后分別計算兩種預(yù)報方法下的實測值和預(yù)測值的偏差并且求對應(yīng)的總偏差平方和，根據(jù)兩個總偏差平方和的大小來得出兩種方法的優(yōu)劣比較，通過Matlab程序的實現(xiàn)得到結(jié)果為：第一種方法比第二種方法好。12/31/202260黃建華制作例5.1雨量預(yù)報方法的評價模型的分析：12/26/2022例5.1雨量預(yù)報方法的評價對于問題二，我們認為公眾的感受主要體現(xiàn)在預(yù)測與實際之間偏差的大小程度，然而，實際測量值未知，因此我們考慮在分級預(yù)報中加入準確概率的方式來實現(xiàn)公眾滿意度的提高。我們認為雨量的實測值應(yīng)該在預(yù)測值點處服從正態(tài)分布，通過合理的假設(shè)和推導(dǎo)，我們得到：正態(tài)分布的均值可以取為雨量預(yù)測值，不同雨量分級區(qū)間上的方差可以近似取為對應(yīng)區(qū)間上的總偏差平方和的平均值。運用Matlab軟件編程，可以實現(xiàn)對每一個預(yù)報數(shù)據(jù)的預(yù)報內(nèi)容的改變，并且通過幾個不同的數(shù)據(jù)體現(xiàn)程序運行所得到的結(jié)果。最后，對模型的缺點進行了討論和改進。

12/31/202261黃建華制作例5.1雨量預(yù)報方法的評價對于問題二，我們認為公眾的感受主例5.1雨量預(yù)報方法的評價模型的假設(shè)假設(shè)測量雨量的工具正常不受任何因素的影響，且所得數(shù)據(jù)真實可靠；測量點所在位置的等距網(wǎng)格點均視為質(zhì)點；符號說明x(i,j)表示第j種方法下的第i個實測點處的雨量實測值y(i,j)表示第j種方法下的第i個實測點處的雨量預(yù)測估計值i=1,2,3…….,91*164,j=1,2e(i,j)表示在第j種方法下的第i個實測點處的偏差Y表示實際值的隨機變量y*表示預(yù)報值ST2(j)表示在第j種方法下的所有實測點處的總偏差平方和12/31/202262黃建華制作例5.1雨量預(yù)報方法的評價模型的假設(shè)12/26/20226例5.1雨量預(yù)報方法的評價問題一模型分析：考察下載的lon.dat、lat.dat、020<日期>.SIX及f<日期><i>_dis<j>等數(shù)據(jù)文件。我們可以看到雨量預(yù)報的網(wǎng)格點有53*47個。而雨量實測點只有91個，且分布不均勻，因此我們不能直接引用數(shù)據(jù)文件進行實測數(shù)據(jù)與預(yù)測數(shù)據(jù)之間偏差的計算。顯然首先要對所提供的數(shù)據(jù)文件進行處理，由于所引用的文件比較多，而且具有一定的規(guī)律性，所以數(shù)據(jù)文件的載入可用MATLAB命令load和循環(huán)結(jié)構(gòu)來實現(xiàn)。當(dāng)有了具體的實測點位置對應(yīng)的實測數(shù)據(jù)yi及預(yù)測數(shù)據(jù)yi*之后,我們可以求出每種方法下所有實測點位置上的總偏差平方和。然后比較兩個總偏差平方和的大小，就可以判斷兩種方法的優(yōu)劣性。

12/31/202263黃建華制作例5.1雨量預(yù)報方法的評價問題一模型分析：12/26/20例5.1雨量預(yù)報方法的評價問題一模型建立：

通過Matlab程序編程計算，具體過程如下：（1）使用MATLAB命令load和循環(huán)結(jié)構(gòu)將lot.dat、lat.dat、f6181_dis1…、f7304_dis2、020618.six…、020730.six等文件輸入MATLAB程序中備用，其中020618.six…、020730.six等文件系統(tǒng)自動在前面添上X。（2）從X020618文件中取第二、第三列作為實測點位置。其中第一列對應(yīng)實測點的緯度，第二列對應(yīng)實測點的經(jīng)度。分別用lat、lon中的數(shù)據(jù)作為網(wǎng)格點的橫坐標和縱坐標。依次取f6181_dis1……,f7304_dis1,f6181_dis2,……,f7304_dis2文件對應(yīng)點數(shù)據(jù)作為豎坐標，張成一個預(yù)報值數(shù)據(jù)曲面，用曲面擬合命令griddata得出對應(yīng)實測點處的預(yù)報估計值，得到一個12/31/202264黃建華制作例5.1雨量預(yù)報方法的評價問題一模型建立：12/26/2例5.1雨量預(yù)報方法的評價91*328維的對應(yīng)實測點處的預(yù)報估計值矩陣yczjz。（3）依次從X020618，…….X020730文件中取第四到第七列組成所有日期，所有時段下的實測數(shù)據(jù)矩陣sczjz。（4）將sczjz分別與yczjz前164列、后164列對應(yīng)元素相減，得出第一種方法下和第二種方法下所有實測點處的偏差，組成矩陣f1pcz和f2pcz。（5）分別求矩陣f1pcz和f2pcz所有元素的平方之和，賦予變量pcpfh(1)和pcpfh(2)。（6）比較變量pcpfh(1)、pcpfh(2)的大小，得出兩種方法優(yōu)劣的比較。各過程的MATLAB程序?qū)崿F(xiàn)，可分別由M-文件sjzr.m、qsj.m、pcbj.m運行得到。

12/31/202265黃建華制作例5.1雨量預(yù)報方法的評價91*328維的對應(yīng)實測點處的預(yù)例5.1雨量預(yù)報方法的評價程序sjzr.m：

rqjz=[618:628,701:730];%產(chǎn)生1*41維日期行向量loadlat.DAT;%載入53*47維緯度矩陣loadlon.DAT;%載入53*47維經(jīng)度矩陣forrq=rqjz%用日期行向量作為循環(huán)向量load(['020',int2str(rq),'.SIX'])%載入53*47維實測值矩陣,系統(tǒng)自動在文件名前加大寫字母Xforsd=1:4%指定時段循環(huán)向量forff=1:2%指定方法循環(huán)向量load(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)])%載入53*47維預(yù)測值矩陣，共328個endendend

12/31/202266黃建華制作例5.1雨量預(yù)報方法的評價程序sjzr.m：12/26/例5.1雨量預(yù)報方法的評價將文件全部拷貝在默認的work文件加下，運行程序sjzr.m：

12/31/202267黃建華制作例5.1雨量預(yù)報方法的評價將文件全部拷貝在默認的work文例5.1雨量預(yù)報方法的評價12/31/202268黃建華制作例5.1雨量預(yù)報方法的評價12/26/202212黃建華制例5.1雨量預(yù)報方法的評價程序huatu.m%畫一個時點預(yù)測曲面圖與實測散點圖比較clf;%預(yù)先清除別的圖象Z=eval(['f',int2str(618),int2str(1),'_dis',int2str(1)]);%取出某時段預(yù)報值作為豎坐標mesh(lat,lon,Z)%張成預(yù)測曲面axis([27,36,117,125,0,0.4])shadingflatxlabel('緯度');ylabel('經(jīng)度');title('圖象比較');holdonX1=X020618(:,2);Y1=X020618(:,3);Z1=X020618(:,4);%取出實測點數(shù)據(jù)plot3(X1,Y1,Z1,'r*')%畫出實測點數(shù)據(jù)散點圖12/31/202269黃建華制作例5.1雨量預(yù)報方法的評價程序huatu.m%畫一例5.1雨量預(yù)報方法的評價12/31/202270黃建華制作例5.1雨量預(yù)報方法的評價12/26/202214黃建華制例5.1雨量預(yù)報方法的評價程序qsj.mycr=1;%指定預(yù)測值矩陣列標yczjz=zeros(91,328);%預(yù)先指定預(yù)測值矩陣為全零矩陣forff=1:2%指定方法循環(huán)變量forrq=rqjz%指定日期循環(huán)變量forsd=1:4%指定時段循環(huán)變量z0=eval(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)]);%依次取出雨量預(yù)測值數(shù)據(jù)yczjz(:,ycr)=griddata(lat,lon,z0,X020618(:,2),X020618(:,3));%在張成的預(yù)測數(shù)據(jù)曲面上擬合對應(yīng)實測點的雨量預(yù)測值，依次放入預(yù)測值矩陣對應(yīng)列ycr=ycr+1;%預(yù)測值矩陣列標向后一列end12/31/202271黃建華制作例5.1雨量預(yù)報方法的評價程序qsj.m12/26/202例5.1雨量預(yù)報方法的評價endendscr=1;%指定實測值矩陣列標sczjz=zeros(91,164);%預(yù)先指定實測值矩陣為全零矩陣forrq=rqjz%指定日期循環(huán)變量sczjz(:,scr:scr+3)=eval(['X020',int2str(rq),'(:,4:7)']);%依次取出雨量實測值數(shù)據(jù)放入實測值矩陣scr=scr+4;%實測值矩陣列標向后4列end

12/31/202272黃建華制作例5.1雨量預(yù)報方法的評價end12/26/2例5.1雨量預(yù)報方法的評價程序pcbj.mff1=yczjz(:,1:164);%取出預(yù)測值矩陣的前164列ff2=yczjz(:,165:328);%取出預(yù)測值矩陣的后164列f1pcz=sczjz-ff1;%計算方法1下的偏差f2pcz=sczjz-ff2;%計算方法2下的偏差pcpfh=[0,0];%預(yù)先產(chǎn)生1*2維偏差平方和零矩陣forpch=1:91%指定偏差行循環(huán)變量forpcr=1:164%指定偏差列循環(huán)變量pcpfh(1)=pcpfh(1)+f1pcz(pch,pcr).^2;%計算方法1偏差平方和pcpfh(2)=pcpfh(2)+f2pcz(pch,pcr).^2;%計算方法2偏差平方和12/31/202273黃建華制作例5.1雨量預(yù)報方法的評價程序pcbj.m12/26/20例5.1雨量預(yù)報方法的評價endendpcpfhifpcpfh(1)==pcpfh(2)%比較兩種方法下的偏差平方和display('一樣好')elseifpcpfh(1)<pcpfh(2)display('第一種方法好')elsedisplay('第二種方法好')end首先將數(shù)據(jù)文件減壓縮在matlab程序的work文件夾下，然后依次運行M-文件sjzr.m、qsj.m、pcbj.m，得到如下結(jié)論：偏差平方和：226730243880第一種方法比第二種方法好。12/31/202274黃建華制作例5.1雨量預(yù)報方法的評價end12/26/202例5.1雨量預(yù)報方法的評價問題二模型分析：在對問題二的分析中，我們只對方法一相應(yīng)的數(shù)據(jù)進行分析。在現(xiàn)實生活中，公眾對天氣預(yù)報的感受主要來自于對雨量預(yù)報等級和雨量實際等級之間的偏差所帶來的不滿，而不是對雨量等級本身的不滿，因此我們認為，考慮公眾的感受就是要盡量減少所預(yù)報等級和實際等級之間偏差的程度。但在預(yù)報過程中，我們不可能引用來自于還未知的實際數(shù)據(jù)，我們只能從評價方法中提供的預(yù)報數(shù)據(jù)及以前的歷史數(shù)據(jù)來提高預(yù)報的準確性及公眾的感受度。為此我們考慮在分級預(yù)報中引入概率。比如：預(yù)報“中雨”時，我們可以采取預(yù)報方式為“中雨，概率為：55%”等。12/31/202275黃建華制作例5.1雨量預(yù)報方法的評價問題二模型分析：12/26/20例5.1雨量預(yù)報方法的評價一般的，雨量實際值是在預(yù)報值附近振動的，而且越接近預(yù)報值的實際值的概率越高，離預(yù)報值越遠的實際值的概率越低，因此我們可以認為：實際值服從以預(yù)報值為中心的正態(tài)分布。假設(shè)實際值隨機變量為Y，預(yù)報值y*則Y~N(y*,σ2)，其中σ2待定?，F(xiàn)假設(shè)對實際值Y進行n次檢測得到樣本為Yi，I=1,2,…..,n,對應(yīng)樣本值為yi，顯然，由概率統(tǒng)計知識可知Yi~N(y*,σ2),且是y*的無偏估計，是σ2的無偏估計。即可用代替y*，s2代替σ2。12/31/202276黃建華制作例5.1雨量預(yù)報方法的評價一般的，雨量實際值是在預(yù)報例5.1雨量預(yù)報方法的評價則：其中ei為每次預(yù)報的偏差。在此，我們顯然不能對同一位置進行n次偏差檢測。但是可以用在同一分級下的實測值與預(yù)測值之間偏差近似代替ei并求σ2值。然后再用各分級下的σ2值求給定預(yù)報值所在對應(yīng)分級下的預(yù)報準確概率P。

12/31/202277黃建華制作例5.1雨量預(yù)報方法的評價則：12/26/202221例5.1雨量預(yù)報方法的評價問題二模型建立：

由上述分析可知，若將無雨、小雨、中雨、大雨、暴雨、大暴雨、特大暴雨等7個等級所對應(yīng)的降雨量區(qū)間記為Ij=[aj,bj],則當(dāng)y*∈Ij時，在第j分級下的均方差為σ2j，Y~N(y*,σ2j),而且Y的實際值y∈Ij的概率為：P=預(yù)報方式為：等級，概率為：P12/31/202278黃建華制作例5.1雨量預(yù)報方法的評價問題二模型建立：12/26/2例5.1雨量預(yù)報方法的評價程序pcbj.mff1=yczjz(:,1:164);%取出預(yù)測值矩陣的前164列ff2=yczjz(:,165:328);%取出預(yù)測值矩陣的后164列f1pcz=sczjz-ff1;%計算方法1下的偏差f2pcz