- 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ì)