官术网_书友最值得收藏!

  • 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)|axb,c(x)≤yd(x)}

需要求解的積分表達(dá)式為:

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

在以上表達(dá)式中wmvn表示的是權(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ù)分析。

主站蜘蛛池模板: 鄂州市| 洛隆县| 靖边县| 新乡县| 寿宁县| 隆子县| 霍城县| 滨州市| 绥芬河市| 广安市| 长汀县| 临西县| 峡江县| 固始县| 左权县| 巧家县| 闻喜县| 郸城县| 驻马店市| 娱乐| 泗阳县| 清新县| 保定市| 吉木萨尔县| 宁夏| 阿鲁科尔沁旗| 扎囊县| 临安市| 庄浪县| 疏附县| 台南市| 娄烦县| 延边| 临猗县| 武隆县| 利津县| 崇明县| 桂平市| 遂昌县| 哈巴河县| 营山县|