第2部分 數據分析
第5章 函數分析和數值運算
本章包括
◆ 函數的零點分析
◆ 一元函數的數值積分
◆ 多重數值積分
◆ 概率分布
◆ 假設檢驗
MATLAB除了可以對矩陣進行分析之外，還可以對函數進行處理。對于常見的高等數學問題，MATLAB都可以有效地進行分析。因此，MATLAB被廣泛應用于數據實驗中。對于微積分、概率論等常見問題，MATLAB提供了對應的命令。讀者可以十分便利地計算和處理對應的問題。
5.1 函數的零點
對于某任意函數f(x)，在求解范圍之內可能有零點，也可能沒有零點，可能只有一個零點也可能有多個甚至是無數個零點。因此，這就給程序求解函數零點增加了很大的難度，沒有可以求解所有函數零點的通用求解命令。在本節中，將簡單討論一元函數和多元函數的零點求解問題。
5.1.1 一元函數的零點
在所有函數中，一元函數是最簡單的，同時也是可以使用MATLAB提供的圖形繪制命令來實現可視化的。因此，在本節中將首先討論一元函數零點的求取方法。在MATLAB中，求解一元函數零點的命令是fzero，其調用格式如下：
◆ x = fzero(fun,x0) 參數fun表示的是一元函數，x0表示求解的初始數值。
◆ [x,fvaI,exitfIag,output] = fzero(fun,x0,options) 參數options的含義是指優化迭代所采用的參數選項，該參數和后面章節中講解到的fsoIve、fminbnd、fminsearch等命令的options都是相同的“模塊”；在輸出參數中，fvaI表示對應的函數值，exitfIag表示的是程序退出的類型，output則是反映優化信息的變量。
例5.1 求函數f(x)=x2sin x-x+1在數值區間[-3,4]中的零點。
step 1 繪制函數的圖形。在MATLAB的命令窗口中輸入以下命令：
%計算函數數值
>> x=[-3:0.1:4];
>>y=sin(x).＊x.^2-x+1;
%繪制函數圖形
>>plot(x,y,'r','LineWidth',1.5)
>>hold on
%添加水平線
>>h=line([-3,4],[0,0]);
%設置直線的寬度和顏色
>>set(h,'LineWidth',1.5)
>>set(h,'color','k')
%設置坐標軸刻度
>>set(gca,'Xtick',[-3:0.5:4])
%添加圖形標題和坐標軸名稱
>>title('The zero of function')
>> grid
>> xlabel('x')
>> ylabel('f(x)')
step 2 查看圖形。輸入以上程序代碼后，按“Enter”鍵，得到的圖形如圖5.1所示。
圖5.1 函數的圖形
說明
之所以在求解函數零點之前，需要繪制函數的圖形，是為了能夠在后面的步驟中使用fzero命令時，更好地選擇初始數值x0。
step 3 求解函數的零點。在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 查看求解的結果。在命令窗口中輸入計算的變量名稱，得到的結果如下：
x =
-2.5708 -1.6194 2.9142
f =
1.0e-015 ＊
-0.8882 0.2220 -0.8882
從以上結果可以看出，函數 f(x)=x2sin x-x+1在[-3,4]范圍內的三個零點數值解為-2.5708、-1.6194和2.9142。
提示
對于一元多項式函數，MATLAB提供了roots命令來求解多項式函數的所有零點，該命令的基本原理是求解多項式伴隨矩陣的特征值來求解。
5.1.2 多元函數的零點
一般來講，多元函數的零點問題比一元函數的零點問題更難解決，但是當零點大致位置和性質比較好預測時，也可以使用數值方法來搜索到精確的零點。
在MATLAB中，求解多元函數的命令是fsoIve，其具體的調用格式如下：
◆ x = fsoIve(fun,x0) 解非線性方程組的數值解。
◆ [x,fvaI,exitfIag,output] = fsoIve(fun,x0,options) 完整格式。
例5.2 求二元方程組
的零點。
step 1 繪制函數的圖形。在MATLAB的命令窗口中輸入以下命令：
%創建三維圖形的數據網格
x=[-5:0.1:5];
y=x;
[X,Y]=meshgrid(x,y);
%計算三維函數的數值
Z=2＊X-Y-exp(-1＊X);
%繪制曲面圖
surf(X,Y,Z)
%設置照明屬性
shading interp
%添加水平的顏色條
colorbar horiz
%設置圖形的坐標軸刻度屬性
set(gca,'Ztick',[-180:20:20])
set(gca,'ZLim',[-170 20])
%設置透明屬性
alphamap('rampdown')
colormap hot
%添加圖形標題和坐標軸名稱
title('The figure of the function')
xlabel('x')
ylabel('y')
zlabel('z')
step 2 查看圖形。在輸入以上代碼后，按“Enter”鍵，得到函數圖形如圖5.2所示。
圖5.2 函數圖形
step 3 編寫求解函數的M文件。選擇命令窗口菜單欄中的“FiIe”→“New”→“M-FiIe”命令，打開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 求解二元函數的零點。在MATLAB的命令窗口中輸入以下代碼：
x0 = [-5; -5];
options=optimset('Display','iter');
[x,fval] = fsolve(@fsolvefun,x0,options)
step 5 查看求解結果。在輸入以上代碼后，按“Enter”鍵，得到以下結果：
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
說明
從以上結果可以看出，由于原來的二元函數是對稱的，因此所求解的未知數結果是相等的，由于在以上實例中設置了顯示迭代，因此在以上結果中顯示各優化信息。
5.2 數值積分
微積分是高等數學的重要知識，在工程實踐中，微積分有著十分廣泛的應用，因此如何通過計算機實現微積分是十分重要的內容。在MATLAB中，可以使用多種方法來實現微積分的運算：數值積分、符號積分、樣條積分和SimuIink模擬積分等。在本章中，主要介紹數值積分和樣條積分，并輔以介紹符號積分和SimuIink積分等方法。
5.2.1 一元函數的數值積分
在MATLAB中，對一元函數進行數值積分的命令是quad和quadI。一般來講，quadI命令比quad命令更加有效，它們的主要功能在于計算閉型數值積分，其對應的詳細調用格式如下：
◆ q = quad(fun,a,b,toI,trace) 采用遞推自適應Simpson法計算積分。
◆ q = quadI(fun,a,b,toI,trace) 采用遞推自適應Lobatto法計算積分。
下面詳細介紹這兩個函數的參數含義：
◆ fun：被積函數，可以是字符串、內聯函數、M函數文件名稱的函數句柄。被積函數中一般使用x作為自變量。
◆ a、b：被積函數的上限和下限，必須都是確定的數值。
◆ toI：標量，控制絕對誤差，默認的數值是精度10-6。
◆ trace：如果該輸入變量的數值不是零，則隨積分的進程逐點繪制被積函數。
說明
在更為完整的調用命令q = quadl(fun,a,b,tol,trace,p1,p2…)中，p1和p2表示的是通過程序向被積函數傳遞的參數。
例5.3 求解積分
的數值。
step 1 分析參數方程。根據微積分的基礎知識，該積分的數值實際上是某曲線的長度，該函數對應的參數方程如下：
step 2 繪制函數圖形。因此為了了解該積分的數學含義，可以繪制該函數的圖形。在MATLAB的命令窗口中輸入以下程序代碼：
%繪制函數圖形
t = 0:0.1:3＊pi;
h=plot3(sin(2＊t),cos(t),t,'r');
set(h,'LineWidth',1.5)
grid on
step 3 查看函數圖形。在輸入以上程序代碼后，按“Enter”鍵，得到的函數圖形如圖5.3所示。
圖5.3 函數圖形
step 4 編寫被積函數的M文件。選擇命令窗口菜單欄中的“FiIe”→“New”→“M-FiIe”命令，打開M文件編輯器，在其中輸入以下程序代碼：
%編寫被積函數的M文件
function f = hcurve(t)
f = sqrt(4＊cos(2＊t).^2 + sin(t).^2 + 1);
在輸入以上程序代碼后，將代碼保存為“hcurve.m”文件。
step 5 使用quad和quadI命令來求解數值積分。在命令窗口中輸入以下程序代碼：
>> len1=quad(@hcurve,0,3＊pi);
>> len2=quadl(@hcurve,0,3＊pi);
step 6 查看求解結果。在命令窗口中輸入變量名稱，得到求解的結果如下：
len1 =
17.2220
len2 =
17.2220
step 7 設置積分的求解屬性，重新求解數值積分。在命令窗口中輸入以下程序代碼：
>>len3=quad(@hcurve,0,3＊pi,1.e-3,1);
>>len4=quadl(@hcurve,0,3＊pi,1.e-3,1);
step 8 查看求解結果。在輸入以上代碼后，按“Enter”鍵，得到以下結果：
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
…………………………………………………………………………//限于篇幅，省略了部分數據
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
………………………………………………………………………… //限于篇幅，省略了部分數據
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
說明
從以上結果可以看出，當設置了比較小的誤差后，使用quad和quadl命令求解得到的結果精度有很大的差別。
5.2.2 使用SimuIink求解數值積分
例5.4 使用SimuIink求解積分
的數值。
step 1 選擇MATLAB菜單欄中的“FiIe” →“New”→“ModeI”命令，打開模型編輯器，向其中添加對應的模型塊，如圖5.4所示。
圖5.4 添加系統的模塊
step 2 設置系統模塊的屬性。雙擊圖5.4中的“MATLAB Fcn”模塊，打開對應的模塊屬性對話框，在其中設置被積函數表達式，如圖5.5所示。
step 3 設置系統仿真的時間，運行系統仿真。將系統仿真的時間設置為3π，然后運行仿真，得到的結果如圖5.6所示。
step 4 查看被積函數的圖形。從圖5.6中可以看出，通過SimuIink積分得到的結果是17.22。同時，可以查看被積函數的圖形，如圖5.7所示。
圖5.5 設置系統模塊的屬性
圖5.6 查看仿真的結果
圖5.7 被積函數的圖形
5.2.3 求解瑕積分
例5.5 求解積分
的數值。
step 1 分析積分的問題。在以上積分表達式中，由于在x=0時，被積函數
的數值趨近于無窮大，也就是
。因此，以上積分屬于瑕積分或者稱為開型積分，這和以上有限積分在性質上有比較大的差異。同時根據基礎的高等數學知識，該積分數值為
。
step 2 求解積分數值。在MATLAB的命令窗口中輸入以下代碼：
>>f=inline('sqrt(log(1./x))','x');
>> lnfq=quad(f,0,1);
>> lnfql=quadl(f,0,1);
step 3 查看求解結果。在命令窗口中輸入變量名稱，得到求解的結果如下：
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 使用符號方法來求解。在MATLAB的命令窗口中輸入以下代碼：
>> f=inline('sqrt(log(1/x))','x');
>> syms x
>> Is=vpa(int('sqrt(log(1/x))','x',0,1))
step 5 查看求解結果。在命令窗口中輸入變量名稱，得到求解的結果如下：
Warning: Explicit integral could not be found.
> In sym.int at 58
In char.int at 9
Is =
.88622692545275801364908374167057
提示
從以上結果可以看出，使用符號方法求解得到的結果和數值積分得到的結果相同，精度都是可以接受的，但是，符號運算比數值運算更加占用系統資源。
5.2.4 矩形區域的多重數值積分
多重數值積分可以認為是一元函數積分的推廣和延伸，但是情況比一元函數要復雜，在本節中，將主要介紹如何在MATLAB中計算二重數值積分。
在MATLAB中，計算二重數值積分的命令為dbIquad，其具體的調用格式如下：
q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method)
這些參數中，fun表示的是積分函數，xmin,xmax表示的是變量x的上下限，ymin，ymax表示的是變量y的上下限；toI表示的是積分絕對誤差，默認數值為10-8；method表示的是積分方法的選項，默認選項為@quad，可以選擇@quadI或者自己定義的積分函數句柄。
例5.6 求解積分
(y sin x+x cos y)dxdy的數值。
step 1 求解積分數值。在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 查看求解結果。輸入以上代碼后，按“Enter”鍵，得到求解的結果如下：
result =
-9.8696
step 3 使用符號運算求解積分數值。
>>syms x y
>> result1=vpa(int(int((y＊sin(x)+x＊cos(y)),x,pi,2＊pi),y,0,pi))
step 4 查看求解結果。輸入以上代碼后，按“Enter”鍵，得到求解的結果如下：
result1 =
-9.8696044010893586188344909998761
5.2.5 變量區域的多重數值積分
前面所介紹的都是固定數值的二重積分運算方法，但是在實際應用中，二重積分并不都是矩形計算區域，在計算區域中會包含變量表達式，也就是說，積分區域可以表示成為以下形式：
R= {(x,y)|a≤x≤b,c(x)≤y≤d(x)}
需要求解的積分表達式為：
對于以上積分表達式，進行數值計算的表達式為：
在以上表達式中wm和vn表示的是權重，取決于一維積分方法。關于二重積分的數值分析的表達式如圖5.8所示。
圖5.8 二重積分數據點
根據圖5.8中數據點區域，需要自行編寫M文件，來計算以上數值積分，在本節中，將使用一個簡單的例子說明如何計算二重數值積分。
例5.7 求解積分
的數值。
step 1 編寫一維數值積分的M文件。選擇命令窗口菜單欄中的“FiIe”→“New”→“M-FiIe”命令，打開M文件編輯器，輸入以下程序代碼：
function INTf=smpsns_fxy(f,x,c,d,N)
%函數f(x,y)的一維數值積分數值;
%對應的積分區域是Ry={c<=y<=d}
%當用戶沒有輸入函數中的N參數時，默認值為100
if nargin<5
N=100;
end
%當參數c=d或者參數N=0，則返回積分數值為0
if abs(d-c)<eps|N<=0
INTf=0;
return;
end
%如果參數N是奇數，則將其加1，變成偶數
if mod(N,2)～=0
N=N+1;
end
%計算單位高度數值
h=(d-c)/N;
%計算節點的y軸坐標值
y=c+[0:N]＊h;
%計算節點的積分函數數值
fxy=feval(f,x,y);
%確定積分的限制范圍
fxy(find(fxy==inf))=realmax;
fxy(find(fxy==-inf))=-realmax;
%計算奇數和偶數節點的x坐標數值
kodd=2:2:N;
keven=3:2:N-1;
%根據積分公式得出積分數值
INTf=h/3＊(fxy(1)+fxy(N+1)+4＊sum(fxy(kodd))+2＊sum(fxy(keven)));
在輸入以上程序代碼后，將代碼保存為“smpsns_fxy.m”文件。
step 2 編寫二重數值積分的M文件。選擇命令窗口菜單欄中的“FiIe”→“New”→“M-FiIe”命令，打開M文件編輯器，輸入以下程序代碼：
function INTfxy=int2s(f,a,b,c,d,M,N)
% 被積函數f(x,y)的二重積分數值
% 積分區域為R={(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;
%判斷參數c是否是數值
%如果c是數值，將積分限制設置為數值c
%如果c不是數值，則將積分顯示設置為函數表達式
if isnumeric(c)
cx(m)=c;
else
cx(m)=feval(c,x(m));
end
%判斷參數d是否是數值
%如果d是數值，將積分限制設置為數值c
%如果d不是數值，則將積分顯示設置為積分表達式
if isnumeric(d)
dx(m)=d;
else
dx(m)=feval(d,x(m));
end
%重復和參數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
%根據Simpson法則計算各個節點的數值
for m=1:M+1
sx(m)=smpsns_fxy(f,x(m),cx(m),dx(m),Nx(m));
end
%計算奇數和偶數的節點
kodd=2:2:M;
keven=3:2:M-1;
%計算積分數值
INTfxy=hx/3＊(sx(1)+sx(M+1)+4＊sum(sx(kodd))+2＊sum(sx(keven)));
在輸入以上程序代碼后，將代碼保存為“int2s.m”文件。
step 3 進行二重積分計算。在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 查看求解結果。在命令窗口中輸入變量名稱，然后按“Enter”鍵，得到求解的結果如下：
>> Vs1
Vs1 =
1.0470
>> Vs2
Vs2 =
1.0470
>> error1
error1 =
-1.5315e-004
>> error2
error2 =
-1.9685e-004
說明
在以上程序結果中，Vs1和Vs2是分別通過不同的計算方法得到的結果，error1和error2是計算結果和真實結果之間的誤差。
step 5 繪制函數圖形。在MATLAB的命令窗口中輸入以下代碼：
>> surf(X,Y,Z)
>>shading interp
>> colormap hsv
>> colorbar horiz
step 6 查看圖形。輸入以上程序代碼后，按“Enter”鍵，得到的結果如圖5.9所示。
圖5.9 函數圖形
5.3 概率論和數理統計
在本節中，將主要介紹在MATLAB中運用概率論和數理統計的方法，主要的內容包括概率分布、數理統計和假設檢驗等。在每個具體的小節中，分別介紹如何在MATLAB中運用相關知識，對于具體的背景知識，請讀者查看對應的書籍。
5.3.1 雙變量的概率分布
概率分布是概率論和數理統計的基礎知識，在MATLAB中，提供了處理常見概率分布的各種命令，包括二項分布、泊松分布、 2χ分布、t 分布等概率分布。這些內容比較簡單，這里就不詳細展開介紹了，感興趣的讀者可以查閱相應的幫助文件。
在本節中將主要介紹如何在MATLAB中處理雙變量或者多變量的概率分布的情況。首先介紹如何處理雙變量t分布（bivariate t distribution）。根據基礎的概率知識，描述雙變量t分布的重要參數是線性相關矩陣κ和自由度η。下面舉例說明如何在MATLAB中顯示多元t分布的圖形。
例5.8 在MATLAB中使用圖形來顯示雙變量t分布（bivariate t distribution），其中兩個變量服從的分布分別為：t(1)和t(5)，也就是說，兩個變量的自由度分別為1和5。下面使用圖形顯示在兩個變量線性相關矩陣κ的不同取值下的分布情況。
step 1 繪制二元概率分布的圖形。在MATLAB的命令窗口中輸入以下代碼：
%設置分布參數
%n代表的是數據點個數
%nu表示的是自由度
%相關系數矩陣為[1 .8; .8 1]
n = 500;
nu = 1;
%產生多元t分布的隨機數值矩陣
T = mvtrnd([1 .8; .8 1], nu, n);
%計算t分布數值的累積概率分布數值
U = tcdf(T,nu);^
%繪制數據點的圖形，并設置圖形的屬性
subplot(2,2,1);
plot(U(:,1),U(:,2),'r.');
title('rho = 0.8');
xlabel('U1');
ylabel('U2')^
%相關系數矩陣為[1 .1; .1 1]
T = mvtrnd([1 .1; .1 1], nu, n);
U = tcdf(T,nu);
%繪制數據點的圖形，并設置圖形的屬性
subplot(2,2,2);
plot(U(:,1),U(:,2),'r.');
title('rho = 0.1');
xlabel('U1');
ylabel('U2');^
%相關系數矩陣為[1-.1; -.1 1]
T = mvtrnd([1-.1; -.1 1], nu, n);
U = tcdf(T,nu);
%繪制數據點的圖形，并設置圖形的屬性
subplot(2,2,3);
plot(U(:,1),U(:,2),'r.');
title('rho = -0.1');
xlabel('U1');
ylabel('U2');^
%相關系數矩陣為[1-.8; -.8 1]
T = mvtrnd([1-.8; -.8 1], nu, n);
U = tcdf(T,nu);
%繪制數據點的圖形，并設置圖形的屬性
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分布的圖形
對于以上程序代碼做如下說明：
◆ 在以上程序代碼中，mvtrnd命令的功能是從多元t分布中產生隨機數據矩陣，關于其具體的用法，可以查看對應的幫助文件。
◆ tcdf命令的功能是產生t分布的累積概率數值，具體的用法請查閱相應的幫助文件。
提示
從圖5.10中可以看出，兩個變量之間的線性相關矩陣系數值的絕對值越接近1，兩個變量的分布越為緊密；而兩個變量之間的線性相關矩陣系數值的絕對值越接近0，兩個變量的分布越為稀疏。
5.3.2 不同概率分布
在MATLAB中，除了繪制兩個相同分布變量之外，還可以繪制兩個不同隨機分布的變量的數據分布圖，下面舉例詳細說明。
例5.9 兩個相關隨機變量，分別服從Gamma分布和t分布，兩個變量相互獨立，且具體的隨機變量參數為Gamma(2,1)和t(5)，在MATLAB中繪制兩個變量的數據分布圖形。
step 1 繪制二元概率分布的圖形。在MATLAB的命令窗口中輸入以下代碼：
subplot(1,1,1);
%設置概率分布的參數
n = 1000;
rho = .7;
nu = 1;
%產生多元t分布的隨機數值矩陣
T = mvtrnd([1 rho; rho 1], nu, n);
%計算t分布數值的累積概率分布數值
U = tcdf(T,nu);
%產生兩個概率分布的數值
X = [gaminv(U(:,1),2,1) tinv(U(:,2),5)];
@@@
%計算兩個直方圖的數值
[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 兩個獨立隨機變量的圖形
5.3.3 數據分布分析
在實際應用中，用戶經常需要根據有限的試驗數據，推測該樣本數據所滿足的數據分布情況，在本節中，將使用一個簡單的實例來演示如何在MATLAB中推測數據分布情況。
例5.10通過命令產生滿足t分布的多元變量，然后使用自定義的概率密度函數來推測兩個變量滿足的多元正態分布N(μ1,μ2,σ1,σ2)，其中參數分別表示均值和標準方差；然后繪制圖形來驗證這種推測是否正確。
step 1 繪制二元概率分布的圖形。在MATLAB的命令窗口中輸入以下代碼：
%產生隨機數據
x = [trnd(20,1,50) trnd(4,1,100)+3];
%設置混合概率密度函數
pdf_normmixture = @(x,p,mu1,mu2,sigma1,sigma2) ...
p＊normpdf(x,mu1,sigma1) + (1-p)＊normpdf(x,mu2,sigma2);
%設置參數的數值
pStart = .5;
muStart = quantile(x,[.25 .75]);
sigmaStart = sqrt(var(x) - .25＊diff(muStart).^2);
start = [pStart muStart sigmaStart sigmaStart];
%設置參數的上下限
lb = [0-Inf -Inf 0 0];
ub = [1 Inf Inf Inf Inf];
%設置求解的屬性
options = statset('MaxIter',300, 'MaxFunEvals',600);
paramEsts = mle(x, 'pdf',pdf_normmixture, 'start',start, ...
'lower',lb, 'upper',ub, 'options',options);
%繪制基礎數據的直方圖
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 數據分布圖形
5.3.4 假設檢驗
假設檢驗是數理統計中的一個重要內容，在MATLAB中可以實現多種常見概率分布的假設檢驗，例如單側檢驗、雙側檢驗、均值檢驗、方差檢驗等多種。在本節中，將以幾個簡單的例子來說明如何在MATLAB中進行假設檢驗。
例5.11 對某試驗數據進行平均值的正態分布單側檢驗，總體的標準差已知，并且假設檢驗的置信水平為5%，假設檢驗的平均值為100。
step 1 進行假設檢驗。在MATLAB的命令窗口中輸入以下代碼：
%設置假設檢驗的參數
mu0 = 100;
sig = 20;
N = 16;
%設置假設檢驗的置信水平
alpha = 0.05;
conf = 1-alpha;
%設置正態分布的截斷點
cutoff = norminv(conf, mu0, sig/sqrt(N));
%產生數據點
x = [linspace(90,cutoff), linspace(cutoff,127)];
y = normpdf(x,mu0,sig/sqrt(N));
%繪制正態分布圖形
h1 = plot(x,y);
xhi = [cutoff, x(x>=cutoff)];
yhi = [0, y(x>=cutoff)];
%繪制假設檢驗的面積圖
patch(xhi,yhi,'b');
%設置圖形的標題和坐標軸名稱
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 單側假設檢驗圖形
提示
在以上程序代碼中，首先輸入了關于該假設檢驗的參數，然后根據正態分布反函數來求取滿足假設檢驗條件的均值數值，最后在正態分布曲線中繪制出假設檢驗的圖形。
step 3 修改假設檢驗的條件，進行假設檢驗。修改假設檢驗條件為均值110，重新進行假設檢驗。在MATLAB的命令窗口中輸入以下代碼：
%修改假設檢驗的均值數值
mu1 = 110;
%計算正態分布數據點
y2 = normpdf(x,mu1,sig/sqrt(N));
%繪制正態分布圖形
h2 = line(x,y2,'Color','r');
%繪制假設檢驗的面積圖
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 修改條件后的假設檢驗
step 5 繪制累積概率密度函數的圖形。為了對比上面兩種不同的假設檢驗條件，可以繪制概率密度函數的圖形。在MATLAB的命令窗口中輸入以下程序代碼：
%計算累積概率密度函數
ynull = normcdf(x,mu0,sig/sqrt(N));
yalt = normcdf(x,mu1,sig/sqrt(N));
%繪制累積密度函數圖形
plot(x,ynull,'b-',x,yalt,'r-');
%計算置信條件水平下的反正態分布數值
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 累積概率密度圖形
說明
圖5.15中顯示了在兩種不同的分布下的假設檢驗情況，實際上，可以在同一個圖形中顯示多個不同分布的假設檢驗情況。
step 7 繪制power函數。power函數的定義是假設檢驗中拒絕檢驗的概率，繪制該函數的圖形也可以查看不同的假設檢驗的情況。在MATLAB的命令窗口中輸入以下程序代碼：
%定義需要的power參數數值
DesiredPower = 0.80;
Nvec = 1:30;
cutoff = mu0 + norminv(conf)＊sig./sqrt(Nvec);
%計算假設檢驗的power數值
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函數的圖形
step 9 使用Monte CarIo模擬檢驗power函數的檢驗結果。在MATLAB的命令窗口中輸入以下程序代碼：
%定義Monte Carlo模擬的參數
%nsamples表示的是樣本數值
%smaplenum是樣本中的數據點
nsamples = 400;
samplenum = 1:nsamples;
N=25;
%創建零值矩陣
h0 = zeros(1,nsamples);
h1 = h0;
%進行右側，已知方差條件下均值假設檢驗
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;
%繪制對應的圖形
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檢驗圖形
5.4 小結
在本章中，依次向讀者介紹了函數零點、數值積分和概率論等內容，這些內容是MATLAB進行數值運算的重要部分。在后面的章節中，將向讀者介紹如何使用MATLAB進行數據分析。