- MATLAB寶典
- 陳杰等編著
- 8字
- 2019-03-01 22:12:33
第2部分 數(shù)據(jù)分析
第5章 函數(shù)分析和數(shù)值運(yùn)算
本章包括
◆ 函數(shù)的零點(diǎn)分析
◆ 一元函數(shù)的數(shù)值積分
◆ 多重?cái)?shù)值積分
◆ 概率分布
◆ 假設(shè)檢驗(yàn)
MATLAB除了可以對(duì)矩陣進(jìn)行分析之外,還可以對(duì)函數(shù)進(jìn)行處理。對(duì)于常見(jiàn)的高等數(shù)學(xué)問(wèn)題,MATLAB都可以有效地進(jìn)行分析。因此,MATLAB被廣泛應(yīng)用于數(shù)據(jù)實(shí)驗(yàn)中。對(duì)于微積分、概率論等常見(jiàn)問(wèn)題,MATLAB提供了對(duì)應(yīng)的命令。讀者可以十分便利地計(jì)算和處理對(duì)應(yīng)的問(wèn)題。
5.1 函數(shù)的零點(diǎn)
對(duì)于某任意函數(shù)f(x),在求解范圍之內(nèi)可能有零點(diǎn),也可能沒(méi)有零點(diǎn),可能只有一個(gè)零點(diǎn)也可能有多個(gè)甚至是無(wú)數(shù)個(gè)零點(diǎn)。因此,這就給程序求解函數(shù)零點(diǎn)增加了很大的難度,沒(méi)有可以求解所有函數(shù)零點(diǎn)的通用求解命令。在本節(jié)中,將簡(jiǎn)單討論一元函數(shù)和多元函數(shù)的零點(diǎn)求解問(wèn)題。
5.1.1 一元函數(shù)的零點(diǎn)
在所有函數(shù)中,一元函數(shù)是最簡(jiǎn)單的,同時(shí)也是可以使用MATLAB提供的圖形繪制命令來(lái)實(shí)現(xiàn)可視化的。因此,在本節(jié)中將首先討論一元函數(shù)零點(diǎn)的求取方法。在MATLAB中,求解一元函數(shù)零點(diǎn)的命令是fzero,其調(diào)用格式如下:
◆ x = fzero(fun,x0) 參數(shù)fun表示的是一元函數(shù),x0表示求解的初始數(shù)值。
◆ [x,fvaI,exitfIag,output] = fzero(fun,x0,options) 參數(shù)options的含義是指優(yōu)化迭代所采用的參數(shù)選項(xiàng),該參數(shù)和后面章節(jié)中講解到的fsoIve、fminbnd、fminsearch等命令的options都是相同的“模塊”;在輸出參數(shù)中,fvaI表示對(duì)應(yīng)的函數(shù)值,exitfIag表示的是程序退出的類型,output則是反映優(yōu)化信息的變量。
例5.1 求函數(shù)f(x)=x2sin x-x+1在數(shù)值區(qū)間[-3,4]中的零點(diǎn)。
step 1 繪制函數(shù)的圖形。在MATLAB的命令窗口中輸入以下命令:
%計(jì)算函數(shù)數(shù)值 >> x=[-3:0.1:4]; >>y=sin(x).*x.^2-x+1; %繪制函數(shù)圖形 >>plot(x,y,'r','LineWidth',1.5) >>hold on %添加水平線 >>h=line([-3,4],[0,0]); %設(shè)置直線的寬度和顏色 >>set(h,'LineWidth',1.5) >>set(h,'color','k') %設(shè)置坐標(biāo)軸刻度 >>set(gca,'Xtick',[-3:0.5:4]) %添加圖形標(biāo)題和坐標(biāo)軸名稱 >>title('The zero of function') >> grid >> xlabel('x') >> ylabel('f(x)')
step 2 查看圖形。輸入以上程序代碼后,按“Enter”鍵,得到的圖形如圖5.1所示。

圖5.1 函數(shù)的圖形
說(shuō)明
之所以在求解函數(shù)零點(diǎn)之前,需要繪制函數(shù)的圖形,是為了能夠在后面的步驟中使用fzero命令時(shí),更好地選擇初始數(shù)值x0。
step 3 求解函數(shù)的零點(diǎn)。在MATLAB的命令窗口中輸入以下命令:
>> [x1,f1,exitflag1]=fzero(f,-2.5); >>[x2,f2,exitflag2]=fzero(f,-1.5); >>[x3,f3,exitflag3]=fzero(f,3); >> x=[x1,x2,x3]; >> f=[f1,f2,f3];
step 4 查看求解的結(jié)果。在命令窗口中輸入計(jì)算的變量名稱,得到的結(jié)果如下:
x = -2.5708 -1.6194 2.9142 f = 1.0e-015 * -0.8882 0.2220 -0.8882
從以上結(jié)果可以看出,函數(shù) f(x)=x2sin x-x+1在[-3,4]范圍內(nèi)的三個(gè)零點(diǎn)數(shù)值解為-2.5708、-1.6194和2.9142。
提示
對(duì)于一元多項(xiàng)式函數(shù),MATLAB提供了roots命令來(lái)求解多項(xiàng)式函數(shù)的所有零點(diǎn),該命令的基本原理是求解多項(xiàng)式伴隨矩陣的特征值來(lái)求解。
5.1.2 多元函數(shù)的零點(diǎn)
一般來(lái)講,多元函數(shù)的零點(diǎn)問(wèn)題比一元函數(shù)的零點(diǎn)問(wèn)題更難解決,但是當(dāng)零點(diǎn)大致位置和性質(zhì)比較好預(yù)測(cè)時(shí),也可以使用數(shù)值方法來(lái)搜索到精確的零點(diǎn)。
在MATLAB中,求解多元函數(shù)的命令是fsoIve,其具體的調(diào)用格式如下:
◆ x = fsoIve(fun,x0) 解非線性方程組的數(shù)值解。
◆ [x,fvaI,exitfIag,output] = fsoIve(fun,x0,options) 完整格式。
例5.2 求二元方程組的零點(diǎn)。
step 1 繪制函數(shù)的圖形。在MATLAB的命令窗口中輸入以下命令:
%創(chuàng)建三維圖形的數(shù)據(jù)網(wǎng)格 x=[-5:0.1:5]; y=x; [X,Y]=meshgrid(x,y); %計(jì)算三維函數(shù)的數(shù)值 Z=2*X-Y-exp(-1*X); %繪制曲面圖 surf(X,Y,Z) %設(shè)置照明屬性 shading interp %添加水平的顏色條 colorbar horiz %設(shè)置圖形的坐標(biāo)軸刻度屬性 set(gca,'Ztick',[-180:20:20]) set(gca,'ZLim',[-170 20]) %設(shè)置透明屬性 alphamap('rampdown') colormap hot %添加圖形標(biāo)題和坐標(biāo)軸名稱 title('The figure of the function') xlabel('x') ylabel('y') zlabel('z')
step 2 查看圖形。在輸入以上代碼后,按“Enter”鍵,得到函數(shù)圖形如圖5.2所示。

圖5.2 函數(shù)圖形
step 3 編寫求解函數(shù)的M文件。選擇命令窗口菜單欄中的“FiIe”→“New”→“M-FiIe”命令,打開(kāi)M文件編輯器,在其中輸入以下程序代碼:
function F = fsolvefun (x) F = [2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))];
在輸入以上程序代碼后,將其保存為“fsoIvefun.m”文件。
step 4 求解二元函數(shù)的零點(diǎn)。在MATLAB的命令窗口中輸入以下代碼:
x0 = [-5; -5]; options=optimset('Display','iter'); [x,fval] = fsolve(@fsolvefun,x0,options)
step 5 查看求解結(jié)果。在輸入以上代碼后,按“Enter”鍵,得到以下結(jié)果:
Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius 0 3 47071.2 2.29e+004 1 1 6 12003.4 1 5.75e+003 1 2 9 3147.02 1 1.47e+003 1 3 12 854.452 1 388 1 4 15 239.527 1 107 1 5 18 67.0412 1 30.8 1 6 21 16.7042 1 9.05 1 7 24 2.42788 1 2.26 1 8 27 0.032658 0.759511 0.206 2.5 9 30 7.03149e-006 0.111927 0.00294 2.5 10 33 3.29525e-013 0.00169132 6.36e-007 2.5 Optimization terminated: first-order optimality is less than options.TolFun. x = 0.5671 0.5671 fval = 1.0e-006 * -0.4059 -0.4059
說(shuō)明
從以上結(jié)果可以看出,由于原來(lái)的二元函數(shù)是對(duì)稱的,因此所求解的未知數(shù)結(jié)果是相等的,由于在以上實(shí)例中設(shè)置了顯示迭代,因此在以上結(jié)果中顯示各優(yōu)化信息。
5.2 數(shù)值積分
微積分是高等數(shù)學(xué)的重要知識(shí),在工程實(shí)踐中,微積分有著十分廣泛的應(yīng)用,因此如何通過(guò)計(jì)算機(jī)實(shí)現(xiàn)微積分是十分重要的內(nèi)容。在MATLAB中,可以使用多種方法來(lái)實(shí)現(xiàn)微積分的運(yùn)算:數(shù)值積分、符號(hào)積分、樣條積分和SimuIink模擬積分等。在本章中,主要介紹數(shù)值積分和樣條積分,并輔以介紹符號(hào)積分和SimuIink積分等方法。
5.2.1 一元函數(shù)的數(shù)值積分
在MATLAB中,對(duì)一元函數(shù)進(jìn)行數(shù)值積分的命令是quad和quadI。一般來(lái)講,quadI命令比quad命令更加有效,它們的主要功能在于計(jì)算閉型數(shù)值積分,其對(duì)應(yīng)的詳細(xì)調(diào)用格式如下:
◆ q = quad(fun,a,b,toI,trace) 采用遞推自適應(yīng)Simpson法計(jì)算積分。
◆ q = quadI(fun,a,b,toI,trace) 采用遞推自適應(yīng)Lobatto法計(jì)算積分。
下面詳細(xì)介紹這兩個(gè)函數(shù)的參數(shù)含義:
◆ fun:被積函數(shù),可以是字符串、內(nèi)聯(lián)函數(shù)、M函數(shù)文件名稱的函數(shù)句柄。被積函數(shù)中一般使用x作為自變量。
◆ a、b:被積函數(shù)的上限和下限,必須都是確定的數(shù)值。
◆ toI:標(biāo)量,控制絕對(duì)誤差,默認(rèn)的數(shù)值是精度10-6。
◆ trace:如果該輸入變量的數(shù)值不是零,則隨積分的進(jìn)程逐點(diǎn)繪制被積函數(shù)。
說(shuō)明
在更為完整的調(diào)用命令q = quadl(fun,a,b,tol,trace,p1,p2…)中,p1和p2表示的是通過(guò)程序向被積函數(shù)傳遞的參數(shù)。
例5.3 求解積分的數(shù)值。
step 1 分析參數(shù)方程。根據(jù)微積分的基礎(chǔ)知識(shí),該積分的數(shù)值實(shí)際上是某曲線的長(zhǎng)度,該函數(shù)對(duì)應(yīng)的參數(shù)方程如下:

step 2 繪制函數(shù)圖形。因此為了了解該積分的數(shù)學(xué)含義,可以繪制該函數(shù)的圖形。在MATLAB的命令窗口中輸入以下程序代碼:
%繪制函數(shù)圖形 t = 0:0.1:3*pi; h=plot3(sin(2*t),cos(t),t,'r'); set(h,'LineWidth',1.5) grid on
step 3 查看函數(shù)圖形。在輸入以上程序代碼后,按“Enter”鍵,得到的函數(shù)圖形如圖5.3所示。

圖5.3 函數(shù)圖形
step 4 編寫被積函數(shù)的M文件。選擇命令窗口菜單欄中的“FiIe”→“New”→“M-FiIe”命令,打開(kāi)M文件編輯器,在其中輸入以下程序代碼:
%編寫被積函數(shù)的M文件 function f = hcurve(t) f = sqrt(4*cos(2*t).^2 + sin(t).^2 + 1);
在輸入以上程序代碼后,將代碼保存為“hcurve.m”文件。
step 5 使用quad和quadI命令來(lái)求解數(shù)值積分。在命令窗口中輸入以下程序代碼:
>> len1=quad(@hcurve,0,3*pi); >> len2=quadl(@hcurve,0,3*pi);
step 6 查看求解結(jié)果。在命令窗口中輸入變量名稱,得到求解的結(jié)果如下:
len1 = 17.2220 len2 = 17.2220
step 7 設(shè)置積分的求解屬性,重新求解數(shù)值積分。在命令窗口中輸入以下程序代碼:
>>len3=quad(@hcurve,0,3*pi,1.e-3,1); >>len4=quadl(@hcurve,0,3*pi,1.e-3,1);
step 8 查看求解結(jié)果。在輸入以上代碼后,按“Enter”鍵,得到以下結(jié)果:
9 0.0000000000 2.55958120e+000 4.5159602105 11 0.0000000000 1.27979060e+000 2.1975964146 13 0.0000000000 6.39895300e-001 1.1908151623 15 0.6398952996 6.39895300e-001 0.9939908843 …………………………………………………………………………//限于篇幅,省略了部分?jǐn)?shù)據(jù) 53 8.1449873615 1.27979060e+000 2.1975964146 55 8.1449873615 6.39895300e-001 0.9939908843 57 8.7848826611 6.39895300e-001 1.1908151623 len3 = 17.2099 18 0.0000000000 4.71238898e+000 15.8755795718 23 0.0000000000 4.32369745e-001 1.4707654142 28 0.0000000000 3.96706633e-002 0.1768555623 33 0.0793413265 7.98333951e-002 0.3426629563 ………………………………………………………………………… //限于篇幅,省略了部分?jǐn)?shù)據(jù) 208 8.6393797974 7.98333951e-002 0.1989870400 213 8.7990465876 9.66808141e-002 0.2889072894 218 8.9924082158 9.66808141e-002 0.3638528604 223 9.1857698440 7.98333951e-002 0.3426629563 228 9.3454366343 3.96706633e-002 0.1768555623 len4 = 17.2220
說(shuō)明
從以上結(jié)果可以看出,當(dāng)設(shè)置了比較小的誤差后,使用quad和quadl命令求解得到的結(jié)果精度有很大的差別。
5.2.2 使用SimuIink求解數(shù)值積分
例5.4 使用SimuIink求解積分的數(shù)值。
step 1 選擇MATLAB菜單欄中的“FiIe” →“New”→“ModeI”命令,打開(kāi)模型編輯器,向其中添加對(duì)應(yīng)的模型塊,如圖5.4所示。

圖5.4 添加系統(tǒng)的模塊
step 2 設(shè)置系統(tǒng)模塊的屬性。雙擊圖5.4中的“MATLAB Fcn”模塊,打開(kāi)對(duì)應(yīng)的模塊屬性對(duì)話框,在其中設(shè)置被積函數(shù)表達(dá)式,如圖5.5所示。
step 3 設(shè)置系統(tǒng)仿真的時(shí)間,運(yùn)行系統(tǒng)仿真。將系統(tǒng)仿真的時(shí)間設(shè)置為3π,然后運(yùn)行仿真,得到的結(jié)果如圖5.6所示。
step 4 查看被積函數(shù)的圖形。從圖5.6中可以看出,通過(guò)SimuIink積分得到的結(jié)果是17.22。同時(shí),可以查看被積函數(shù)的圖形,如圖5.7所示。

圖5.5 設(shè)置系統(tǒng)模塊的屬性

圖5.6 查看仿真的結(jié)果

圖5.7 被積函數(shù)的圖形
5.2.3 求解瑕積分
例5.5 求解積分的數(shù)值。
step 1 分析積分的問(wèn)題。在以上積分表達(dá)式中,由于在x=0時(shí),被積函數(shù)的數(shù)值趨近于無(wú)窮大,也就是
。因此,以上積分屬于瑕積分或者稱為開(kāi)型積分,這和以上有限積分在性質(zhì)上有比較大的差異。同時(shí)根據(jù)基礎(chǔ)的高等數(shù)學(xué)知識(shí),該積分?jǐn)?shù)值為
。
step 2 求解積分?jǐn)?shù)值。在MATLAB的命令窗口中輸入以下代碼:
>>f=inline('sqrt(log(1./x))','x'); >> lnfq=quad(f,0,1); >> lnfql=quadl(f,0,1);
step 3 查看求解結(jié)果。在命令窗口中輸入變量名稱,得到求解的結(jié)果如下:
Warning: Divide by zero. > In inlineeval at 13 In inline.subsref at 25 In quad at 62 lnfq = 0.8862 Warning: Divide by zero. > In inlineeval at 13 In inline.feval at 34 In quadl at 64 lnfql = 0.8862
step 4 使用符號(hào)方法來(lái)求解。在MATLAB的命令窗口中輸入以下代碼:
>> f=inline('sqrt(log(1/x))','x'); >> syms x >> Is=vpa(int('sqrt(log(1/x))','x',0,1))
step 5 查看求解結(jié)果。在命令窗口中輸入變量名稱,得到求解的結(jié)果如下:
Warning: Explicit integral could not be found. > In sym.int at 58 In char.int at 9 Is = .88622692545275801364908374167057
提示
從以上結(jié)果可以看出,使用符號(hào)方法求解得到的結(jié)果和數(shù)值積分得到的結(jié)果相同,精度都是可以接受的,但是,符號(hào)運(yùn)算比數(shù)值運(yùn)算更加占用系統(tǒng)資源。
5.2.4 矩形區(qū)域的多重?cái)?shù)值積分
多重?cái)?shù)值積分可以認(rèn)為是一元函數(shù)積分的推廣和延伸,但是情況比一元函數(shù)要復(fù)雜,在本節(jié)中,將主要介紹如何在MATLAB中計(jì)算二重?cái)?shù)值積分。
在MATLAB中,計(jì)算二重?cái)?shù)值積分的命令為dbIquad,其具體的調(diào)用格式如下:
q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method)
這些參數(shù)中,fun表示的是積分函數(shù),xmin,xmax表示的是變量x的上下限,ymin,ymax表示的是變量y的上下限;toI表示的是積分絕對(duì)誤差,默認(rèn)數(shù)值為10-8;method表示的是積分方法的選項(xiàng),默認(rèn)選項(xiàng)為@quad,可以選擇@quadI或者自己定義的積分函數(shù)句柄。
例5.6 求解積分(y sin x+x cos y)dxdy的數(shù)值。
step 1 求解積分?jǐn)?shù)值。在MATLAB的命令窗口中輸入以下代碼:
>>integrnd=@(x,y) y*sin(x)+x*cos(y); >>xmin = pi; >>xmax = 2*pi; >>ymin = 0; >>ymax = pi; >>result = dblquad(integrnd,xmin,xmax,ymin,ymax)
step 2 查看求解結(jié)果。輸入以上代碼后,按“Enter”鍵,得到求解的結(jié)果如下:
result = -9.8696
step 3 使用符號(hào)運(yùn)算求解積分?jǐn)?shù)值。
>>syms x y >> result1=vpa(int(int((y*sin(x)+x*cos(y)),x,pi,2*pi),y,0,pi))
step 4 查看求解結(jié)果。輸入以上代碼后,按“Enter”鍵,得到求解的結(jié)果如下:
result1 = -9.8696044010893586188344909998761
5.2.5 變量區(qū)域的多重?cái)?shù)值積分
前面所介紹的都是固定數(shù)值的二重積分運(yùn)算方法,但是在實(shí)際應(yīng)用中,二重積分并不都是矩形計(jì)算區(qū)域,在計(jì)算區(qū)域中會(huì)包含變量表達(dá)式,也就是說(shuō),積分區(qū)域可以表示成為以下形式:
R= {(x,y)|a≤x≤b,c(x)≤y≤d(x)}
需要求解的積分表達(dá)式為:

對(duì)于以上積分表達(dá)式,進(jìn)行數(shù)值計(jì)算的表達(dá)式為:

在以上表達(dá)式中wm和vn表示的是權(quán)重,取決于一維積分方法。關(guān)于二重積分的數(shù)值分析的表達(dá)式如圖5.8所示。

圖5.8 二重積分?jǐn)?shù)據(jù)點(diǎn)
根據(jù)圖5.8中數(shù)據(jù)點(diǎn)區(qū)域,需要自行編寫M文件,來(lái)計(jì)算以上數(shù)值積分,在本節(jié)中,將使用一個(gè)簡(jiǎn)單的例子說(shuō)明如何計(jì)算二重?cái)?shù)值積分。
例5.7 求解積分的數(shù)值。
step 1 編寫一維數(shù)值積分的M文件。選擇命令窗口菜單欄中的“FiIe”→“New”→“M-FiIe”命令,打開(kāi)M文件編輯器,輸入以下程序代碼:
function INTf=smpsns_fxy(f,x,c,d,N) %函數(shù)f(x,y)的一維數(shù)值積分?jǐn)?shù)值; %對(duì)應(yīng)的積分區(qū)域是Ry={c<=y<=d} %當(dāng)用戶沒(méi)有輸入函數(shù)中的N參數(shù)時(shí),默認(rèn)值為100 if nargin<5 N=100; end %當(dāng)參數(shù)c=d或者參數(shù)N=0,則返回積分?jǐn)?shù)值為0 if abs(d-c)<eps|N<=0 INTf=0; return; end %如果參數(shù)N是奇數(shù),則將其加1,變成偶數(shù) if mod(N,2)~=0 N=N+1; end %計(jì)算單位高度數(shù)值 h=(d-c)/N; %計(jì)算節(jié)點(diǎn)的y軸坐標(biāo)值 y=c+[0:N]*h; %計(jì)算節(jié)點(diǎn)的積分函數(shù)數(shù)值 fxy=feval(f,x,y); %確定積分的限制范圍 fxy(find(fxy==inf))=realmax; fxy(find(fxy==-inf))=-realmax; %計(jì)算奇數(shù)和偶數(shù)節(jié)點(diǎn)的x坐標(biāo)數(shù)值 kodd=2:2:N; keven=3:2:N-1; %根據(jù)積分公式得出積分?jǐn)?shù)值 INTf=h/3*(fxy(1)+fxy(N+1)+4*sum(fxy(kodd))+2*sum(fxy(keven)));
在輸入以上程序代碼后,將代碼保存為“smpsns_fxy.m”文件。
step 2 編寫二重?cái)?shù)值積分的M文件。選擇命令窗口菜單欄中的“FiIe”→“New”→“M-FiIe”命令,打開(kāi)M文件編輯器,輸入以下程序代碼:
function INTfxy=int2s(f,a,b,c,d,M,N) % 被積函數(shù)f(x,y)的二重積分?jǐn)?shù)值 % 積分區(qū)域?yàn)镽={(x,y)|a<=x<=b,c(x)<=y<=d(x)} % 使用的積分方法是Simpson法則 if ceil(M)~=floor(M) hx=M; M=ceil((b-a)/hx); end if mod(M,2)~=0 M=M+1; end hx=(b-a)/M; m=1:M+1; x=a+(m-1)*hx; %判斷參數(shù)c是否是數(shù)值 %如果c是數(shù)值,將積分限制設(shè)置為數(shù)值c %如果c不是數(shù)值,則將積分顯示設(shè)置為函數(shù)表達(dá)式 if isnumeric(c) cx(m)=c; else cx(m)=feval(c,x(m)); end %判斷參數(shù)d是否是數(shù)值 %如果d是數(shù)值,將積分限制設(shè)置為數(shù)值c %如果d不是數(shù)值,則將積分顯示設(shè)置為積分表達(dá)式 if isnumeric(d) dx(m)=d; else dx(m)=feval(d,x(m)); end %重復(fù)和參數(shù)M類似的操作 if ceil(N)~=floor(N) hy=N; Nx(m)=ceil((dx(m)-cx(m))/hy); ind=find(mod(Nx(m),2)~=0); Nx(ind)=Nx(ind)+1; else if mod(N,2)~=0 N=N+1; end Nx(m)=N; end %根據(jù)Simpson法則計(jì)算各個(gè)節(jié)點(diǎn)的數(shù)值 for m=1:M+1 sx(m)=smpsns_fxy(f,x(m),cx(m),dx(m),Nx(m)); end %計(jì)算奇數(shù)和偶數(shù)的節(jié)點(diǎn) kodd=2:2:M; keven=3:2:M-1; %計(jì)算積分?jǐn)?shù)值 INTfxy=hx/3*(sx(1)+sx(M+1)+4*sum(sx(kodd))+2*sum(sx(keven)));
在輸入以上程序代碼后,將代碼保存為“int2s.m”文件。
step 3 進(jìn)行二重積分計(jì)算。在MATLAB的命令窗口中輸入以下代碼:
>> x=[-1:0.05:1]; >> y=[0:0.05:1]; >>[X,Y]=meshgrid(x,y); >>f510=inline('sqrt(max(1-x.*x-y.*y,0))','x','y'); >> Z=f510(X,Y); >> d=inline('sqrt(max(1-x.*x,0))','x'); >> b=1; >> a=-1; >> c=0; >> Vs1=int2s(f510,a,b,c,d,100,100); >> error1=Vs1-pi/3; >> Vs2=int2s(f510,a,b,c,d,0.01,0.01); >> error2=Vs2-pi/3;
step 4 查看求解結(jié)果。在命令窗口中輸入變量名稱,然后按“Enter”鍵,得到求解的結(jié)果如下:
>> Vs1 Vs1 = 1.0470 >> Vs2 Vs2 = 1.0470 >> error1 error1 = -1.5315e-004 >> error2 error2 = -1.9685e-004
說(shuō)明
在以上程序結(jié)果中,Vs1和Vs2是分別通過(guò)不同的計(jì)算方法得到的結(jié)果,error1和error2是計(jì)算結(jié)果和真實(shí)結(jié)果之間的誤差。
step 5 繪制函數(shù)圖形。在MATLAB的命令窗口中輸入以下代碼:
>> surf(X,Y,Z) >>shading interp >> colormap hsv >> colorbar horiz
step 6 查看圖形。輸入以上程序代碼后,按“Enter”鍵,得到的結(jié)果如圖5.9所示。

圖5.9 函數(shù)圖形
5.3 概率論和數(shù)理統(tǒng)計(jì)
在本節(jié)中,將主要介紹在MATLAB中運(yùn)用概率論和數(shù)理統(tǒng)計(jì)的方法,主要的內(nèi)容包括概率分布、數(shù)理統(tǒng)計(jì)和假設(shè)檢驗(yàn)等。在每個(gè)具體的小節(jié)中,分別介紹如何在MATLAB中運(yùn)用相關(guān)知識(shí),對(duì)于具體的背景知識(shí),請(qǐng)讀者查看對(duì)應(yīng)的書籍。
5.3.1 雙變量的概率分布
概率分布是概率論和數(shù)理統(tǒng)計(jì)的基礎(chǔ)知識(shí),在MATLAB中,提供了處理常見(jiàn)概率分布的各種命令,包括二項(xiàng)分布、泊松分布、 2χ分布、t 分布等概率分布。這些內(nèi)容比較簡(jiǎn)單,這里就不詳細(xì)展開(kāi)介紹了,感興趣的讀者可以查閱相應(yīng)的幫助文件。
在本節(jié)中將主要介紹如何在MATLAB中處理雙變量或者多變量的概率分布的情況。首先介紹如何處理雙變量t分布(bivariate t distribution)。根據(jù)基礎(chǔ)的概率知識(shí),描述雙變量t分布的重要參數(shù)是線性相關(guān)矩陣κ和自由度η。下面舉例說(shuō)明如何在MATLAB中顯示多元t分布的圖形。
例5.8 在MATLAB中使用圖形來(lái)顯示雙變量t分布(bivariate t distribution),其中兩個(gè)變量服從的分布分別為:t(1)和t(5),也就是說(shuō),兩個(gè)變量的自由度分別為1和5。下面使用圖形顯示在兩個(gè)變量線性相關(guān)矩陣κ的不同取值下的分布情況。
step 1 繪制二元概率分布的圖形。在MATLAB的命令窗口中輸入以下代碼:
%設(shè)置分布參數(shù) %n代表的是數(shù)據(jù)點(diǎn)個(gè)數(shù) %nu表示的是自由度 %相關(guān)系數(shù)矩陣為[1 .8; .8 1] n = 500; nu = 1; %產(chǎn)生多元t分布的隨機(jī)數(shù)值矩陣 T = mvtrnd([1 .8; .8 1], nu, n); %計(jì)算t分布數(shù)值的累積概率分布數(shù)值 U = tcdf(T,nu);^ %繪制數(shù)據(jù)點(diǎn)的圖形,并設(shè)置圖形的屬性 subplot(2,2,1); plot(U(:,1),U(:,2),'r.'); title('rho = 0.8'); xlabel('U1'); ylabel('U2')^ %相關(guān)系數(shù)矩陣為[1 .1; .1 1] T = mvtrnd([1 .1; .1 1], nu, n); U = tcdf(T,nu); %繪制數(shù)據(jù)點(diǎn)的圖形,并設(shè)置圖形的屬性 subplot(2,2,2); plot(U(:,1),U(:,2),'r.'); title('rho = 0.1'); xlabel('U1'); ylabel('U2');^ %相關(guān)系數(shù)矩陣為[1-.1; -.1 1] T = mvtrnd([1-.1; -.1 1], nu, n); U = tcdf(T,nu); %繪制數(shù)據(jù)點(diǎn)的圖形,并設(shè)置圖形的屬性 subplot(2,2,3); plot(U(:,1),U(:,2),'r.'); title('rho = -0.1'); xlabel('U1'); ylabel('U2');^ %相關(guān)系數(shù)矩陣為[1-.8; -.8 1] T = mvtrnd([1-.8; -.8 1], nu, n); U = tcdf(T,nu); %繪制數(shù)據(jù)點(diǎn)的圖形,并設(shè)置圖形的屬性 subplot(2,2,4); plot(U(:,1),U(:,2),'r.'); title('rho = -0.8'); xlabel('U1'); ylabel('U2')
step 2 查看圖形。在輸入程序代碼后,按“Enter”鍵,得到的圖形如圖5.10所示。

圖5.10 雙變量t分布的圖形
對(duì)于以上程序代碼做如下說(shuō)明:
◆ 在以上程序代碼中,mvtrnd命令的功能是從多元t分布中產(chǎn)生隨機(jī)數(shù)據(jù)矩陣,關(guān)于其具體的用法,可以查看對(duì)應(yīng)的幫助文件。
◆ tcdf命令的功能是產(chǎn)生t分布的累積概率數(shù)值,具體的用法請(qǐng)查閱相應(yīng)的幫助文件。
提示
從圖5.10中可以看出,兩個(gè)變量之間的線性相關(guān)矩陣系數(shù)值的絕對(duì)值越接近1,兩個(gè)變量的分布越為緊密;而兩個(gè)變量之間的線性相關(guān)矩陣系數(shù)值的絕對(duì)值越接近0,兩個(gè)變量的分布越為稀疏。
5.3.2 不同概率分布
在MATLAB中,除了繪制兩個(gè)相同分布變量之外,還可以繪制兩個(gè)不同隨機(jī)分布的變量的數(shù)據(jù)分布圖,下面舉例詳細(xì)說(shuō)明。
例5.9 兩個(gè)相關(guān)隨機(jī)變量,分別服從Gamma分布和t分布,兩個(gè)變量相互獨(dú)立,且具體的隨機(jī)變量參數(shù)為Gamma(2,1)和t(5),在MATLAB中繪制兩個(gè)變量的數(shù)據(jù)分布圖形。
step 1 繪制二元概率分布的圖形。在MATLAB的命令窗口中輸入以下代碼:
subplot(1,1,1); %設(shè)置概率分布的參數(shù) n = 1000; rho = .7; nu = 1; %產(chǎn)生多元t分布的隨機(jī)數(shù)值矩陣 T = mvtrnd([1 rho; rho 1], nu, n); %計(jì)算t分布數(shù)值的累積概率分布數(shù)值 U = tcdf(T,nu); %產(chǎn)生兩個(gè)概率分布的數(shù)值 X = [gaminv(U(:,1),2,1) tinv(U(:,2),5)];
@@@
%計(jì)算兩個(gè)直方圖的數(shù)值 [n1,ctr1] = hist(X(:,1),20); [n2,ctr2] = hist(X(:,2),20); %繪制圖形 subplot(2,2,2); plot(X(:,1),X(:,2),'.'); axis([0 15-10 10]); h1 = gca; title('1000 Simulated Dependent t and Gamma Values'); xlabel('X1~ Gamma(2,1)'); ylabel('X2~ t(5)'); subplot(2,2,4); bar(ctr1,-n1,1); axis([0 15-max(n1)*1.1 0]); axis('off'); h2 = gca; subplot(2,2,1); barh(ctr2,-n2,1); axis([-max(n2)*1.1 0-10 10]); axis('off'); h3 = gca; set(h1,'Position',[0.35 0.35 0.55 0.55]); set(h2,'Position',[.35 .1 .55 .15]); set(h3,'Position',[.1 .35 .15 .55]); colormap([.8 .8 1]);
step 2 查看圖形。在輸入程序代碼后,按“Enter”鍵,得到的圖形如圖5.11所示。

圖5.11 兩個(gè)獨(dú)立隨機(jī)變量的圖形
5.3.3 數(shù)據(jù)分布分析
在實(shí)際應(yīng)用中,用戶經(jīng)常需要根據(jù)有限的試驗(yàn)數(shù)據(jù),推測(cè)該樣本數(shù)據(jù)所滿足的數(shù)據(jù)分布情況,在本節(jié)中,將使用一個(gè)簡(jiǎn)單的實(shí)例來(lái)演示如何在MATLAB中推測(cè)數(shù)據(jù)分布情況。
例5.10通過(guò)命令產(chǎn)生滿足t分布的多元變量,然后使用自定義的概率密度函數(shù)來(lái)推測(cè)兩個(gè)變量滿足的多元正態(tài)分布N(μ1,μ2,σ1,σ2),其中參數(shù)分別表示均值和標(biāo)準(zhǔn)方差;然后繪制圖形來(lái)驗(yàn)證這種推測(cè)是否正確。
step 1 繪制二元概率分布的圖形。在MATLAB的命令窗口中輸入以下代碼:
%產(chǎn)生隨機(jī)數(shù)據(jù) x = [trnd(20,1,50) trnd(4,1,100)+3]; %設(shè)置混合概率密度函數(shù) pdf_normmixture = @(x,p,mu1,mu2,sigma1,sigma2) ... p*normpdf(x,mu1,sigma1) + (1-p)*normpdf(x,mu2,sigma2); %設(shè)置參數(shù)的數(shù)值 pStart = .5; muStart = quantile(x,[.25 .75]); sigmaStart = sqrt(var(x) - .25*diff(muStart).^2); start = [pStart muStart sigmaStart sigmaStart]; %設(shè)置參數(shù)的上下限 lb = [0-Inf -Inf 0 0]; ub = [1 Inf Inf Inf Inf]; %設(shè)置求解的屬性 options = statset('MaxIter',300, 'MaxFunEvals',600); paramEsts = mle(x, 'pdf',pdf_normmixture, 'start',start, ... 'lower',lb, 'upper',ub, 'options',options); %繪制基礎(chǔ)數(shù)據(jù)的直方圖 bins = -2.5:.5:7.5; h = bar(bins,histc(x,bins)/(length(x)*.5),'histc'); set(h,'FaceColor',[.9 .9 .9]); xgrid = linspace(1.1*min(x),1.1*max(x),200); %繪制概率密度圖形 pdfgrid = pdf_normmixture(xgrid,paramEsts(1),paramEsts(2),paramEsts(3), paramEsts(4),paramEsts(5)); hold on; plot(xgrid,pdfgrid,'-'); hold off xlabel('x'); ylabel('Probability Density');
step 2 查看圖形。在輸入程序代碼后,按“Enter”鍵,得到的圖形如圖5.12所示。

圖5.12 數(shù)據(jù)分布圖形
5.3.4 假設(shè)檢驗(yàn)
假設(shè)檢驗(yàn)是數(shù)理統(tǒng)計(jì)中的一個(gè)重要內(nèi)容,在MATLAB中可以實(shí)現(xiàn)多種常見(jiàn)概率分布的假設(shè)檢驗(yàn),例如單側(cè)檢驗(yàn)、雙側(cè)檢驗(yàn)、均值檢驗(yàn)、方差檢驗(yàn)等多種。在本節(jié)中,將以幾個(gè)簡(jiǎn)單的例子來(lái)說(shuō)明如何在MATLAB中進(jìn)行假設(shè)檢驗(yàn)。
例5.11 對(duì)某試驗(yàn)數(shù)據(jù)進(jìn)行平均值的正態(tài)分布單側(cè)檢驗(yàn),總體的標(biāo)準(zhǔn)差已知,并且假設(shè)檢驗(yàn)的置信水平為5%,假設(shè)檢驗(yàn)的平均值為100。
step 1 進(jìn)行假設(shè)檢驗(yàn)。在MATLAB的命令窗口中輸入以下代碼:
%設(shè)置假設(shè)檢驗(yàn)的參數(shù) mu0 = 100; sig = 20; N = 16; %設(shè)置假設(shè)檢驗(yàn)的置信水平 alpha = 0.05; conf = 1-alpha; %設(shè)置正態(tài)分布的截?cái)帱c(diǎn) cutoff = norminv(conf, mu0, sig/sqrt(N)); %產(chǎn)生數(shù)據(jù)點(diǎn) x = [linspace(90,cutoff), linspace(cutoff,127)]; y = normpdf(x,mu0,sig/sqrt(N)); %繪制正態(tài)分布圖形 h1 = plot(x,y); xhi = [cutoff, x(x>=cutoff)]; yhi = [0, y(x>=cutoff)]; %繪制假設(shè)檢驗(yàn)的面積圖 patch(xhi,yhi,'b'); %設(shè)置圖形的標(biāo)題和坐標(biāo)軸名稱 title('Distribution of sample mean, N=16'); xlabel('Sample mean'); ylabel('Density'); text(96,.01,sprintf('Reject if mean>%.4g\nProb = 0.05',cutoff),'Color','b')
step 2 查看圖形。在輸入程序代碼后,按“Enter”鍵,得到的圖形如圖5.13所示。

圖5.13 單側(cè)假設(shè)檢驗(yàn)圖形
提示
在以上程序代碼中,首先輸入了關(guān)于該假設(shè)檢驗(yàn)的參數(shù),然后根據(jù)正態(tài)分布反函數(shù)來(lái)求取滿足假設(shè)檢驗(yàn)條件的均值數(shù)值,最后在正態(tài)分布曲線中繪制出假設(shè)檢驗(yàn)的圖形。
step 3 修改假設(shè)檢驗(yàn)的條件,進(jìn)行假設(shè)檢驗(yàn)。修改假設(shè)檢驗(yàn)條件為均值110,重新進(jìn)行假設(shè)檢驗(yàn)。在MATLAB的命令窗口中輸入以下代碼:
%修改假設(shè)檢驗(yàn)的均值數(shù)值 mu1 = 110; %計(jì)算正態(tài)分布數(shù)據(jù)點(diǎn) y2 = normpdf(x,mu1,sig/sqrt(N)); %繪制正態(tài)分布圖形 h2 = line(x,y2,'Color','r'); %繪制假設(shè)檢驗(yàn)的面積圖 yhi = [0, y2(x>=cutoff)]; patch(xhi,yhi,'r','FaceAlpha',0.25); %添加圖形的提示信息 P = 1- normcdf(cutoff,mu1,sig/sqrt(N)); text(115,.06,sprintf('Reject if T>%.4g\nProb = %.2g',cutoff,P),'Color',[1 0 0]); legend([h1 h2],'Null hypothesis','Alternative hypothesis');
step 4 查看圖形。在輸入程序代碼后,按“Enter”鍵,得到的圖形如圖5.14所示。

圖5.14 修改條件后的假設(shè)檢驗(yàn)
step 5 繪制累積概率密度函數(shù)的圖形。為了對(duì)比上面兩種不同的假設(shè)檢驗(yàn)條件,可以繪制概率密度函數(shù)的圖形。在MATLAB的命令窗口中輸入以下程序代碼:
%計(jì)算累積概率密度函數(shù) ynull = normcdf(x,mu0,sig/sqrt(N)); yalt = normcdf(x,mu1,sig/sqrt(N)); %繪制累積密度函數(shù)圖形 plot(x,ynull,'b-',x,yalt,'r-'); %計(jì)算置信條件水平下的反正態(tài)分布數(shù)值 zval = norminv(conf); cutoff = mu0 + zval * sig / sqrt(N); %繪制圖形 line([90,cutoff,cutoff],[conf, conf, 0],'LineStyle',':'); msg = sprintf(' Cutoff = \\mu_0 + %.2g\\sigma / \\surd{n}',zval); text(cutoff,.15,msg,'Color','b'); text(min(x),conf,sprintf(' %g%% test',100*alpha),'Color','b',... 'verticalalignment','top') palt = normcdf(cutoff,mu1,sig/sqrt(N)); line([90,cutoff],[palt,palt],'Color','r','LineStyle',':'); text(91,palt+.02,sprintf(' Power is 1-%.2g = %.2g',palt,1-palt),'Color', [1 0 0]);
step 6 查看圖形。在輸入程序代碼后,按“Enter”鍵,得到的圖形如圖5.15所示。

圖5.15 累積概率密度圖形
說(shuō)明
圖5.15中顯示了在兩種不同的分布下的假設(shè)檢驗(yàn)情況,實(shí)際上,可以在同一個(gè)圖形中顯示多個(gè)不同分布的假設(shè)檢驗(yàn)情況。
step 7 繪制power函數(shù)。power函數(shù)的定義是假設(shè)檢驗(yàn)中拒絕檢驗(yàn)的概率,繪制該函數(shù)的圖形也可以查看不同的假設(shè)檢驗(yàn)的情況。在MATLAB的命令窗口中輸入以下程序代碼:
%定義需要的power參數(shù)數(shù)值 DesiredPower = 0.80; Nvec = 1:30; cutoff = mu0 + norminv(conf)*sig./sqrt(Nvec); %計(jì)算假設(shè)檢驗(yàn)的power數(shù)值 power = 1- normcdf(cutoff, mu1, sig./sqrt(Nvec)); %繪制圖形 plot(Nvec,power,'bo-',[0 30],[DesiredPower DesiredPower],'k:'); xlabel('N = sample size'); ylabel('Power') title('Power function for the alternative: \mu = 110')
step 8 查看圖形。在輸入程序代碼后,按“Enter”鍵,得到的圖形如圖5.16所示。

圖5.16 power函數(shù)的圖形
step 9 使用Monte CarIo模擬檢驗(yàn)power函數(shù)的檢驗(yàn)結(jié)果。在MATLAB的命令窗口中輸入以下程序代碼:
%定義Monte Carlo模擬的參數(shù) %nsamples表示的是樣本數(shù)值 %smaplenum是樣本中的數(shù)據(jù)點(diǎn) nsamples = 400; samplenum = 1:nsamples; N=25; %創(chuàng)建零值矩陣 h0 = zeros(1,nsamples); h1 = h0; %進(jìn)行右側(cè),已知方差條件下均值假設(shè)檢驗(yàn) for j=1:nsamples Z0 = normrnd(mu0,sig,N,1); h0(j) = ztest(Z0,mu0,sig,alpha,'right'); Z1 = normrnd(mu1,sig,N,1); h1(j) = ztest(Z1,mu0,sig,alpha,'right'); end p0 = cumsum(h0) ./ samplenum; p1 = cumsum(h1) ./ samplenum; %繪制對(duì)應(yīng)的圖形 plot(samplenum,p0,'b-',samplenum,p1,'r-') xlabel('Sample number'); ylabel('Proportion significant') title('Verification of power computation') legend('Null hypothesis','Alternative hypothesis','Location','East')
step 10 查看圖形。在輸入程序代碼后,按“Enter”鍵,得到的圖形如圖5.17所示。

圖5.17 Monte CarIo檢驗(yàn)圖形
5.4 小結(jié)
在本章中,依次向讀者介紹了函數(shù)零點(diǎn)、數(shù)值積分和概率論等內(nèi)容,這些內(nèi)容是MATLAB進(jìn)行數(shù)值運(yùn)算的重要部分。在后面的章節(jié)中,將向讀者介紹如何使用MATLAB進(jìn)行數(shù)據(jù)分析。
- Mastering Mesos
- 高效能辦公必修課:Word圖文處理
- 零起步輕松學(xué)單片機(jī)技術(shù)(第2版)
- Word 2003、Excel 2003、PowerPoint 2003上機(jī)指導(dǎo)與練習(xí)
- 腦動(dòng)力:Linux指令速查效率手冊(cè)
- 中文版Photoshop CS5數(shù)碼照片處理完全自學(xué)一本通
- Mastering VMware vSphere 6.5
- Learning Social Media Analytics with R
- 來(lái)吧!帶你玩轉(zhuǎn)Excel VBA
- HTML5 Canvas Cookbook
- 新編計(jì)算機(jī)圖形學(xué)
- Dreamweaver CS6中文版多功能教材
- 格蠹匯編
- 嵌入式Linux系統(tǒng)實(shí)用開(kāi)發(fā)
- 重估:人工智能與賦能社會(huì)