第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进行数据分析。