概率论与数理统计是研究和应用随机现象统计规律性的一门数学科学.其应用十分广泛,几乎遍及所有科学领域、工农业生产和国民经济各部门。本章将利用Matlab来解决概率统计学中的概率分布、数字特征、参数估计以及假设检验等问题.
8.1 数据分析
8。1。1 几种均值
在给定的一组数据中,要进行各种均值的计算，在Matlab中可由以下函数实现。
mean 算术平均值函数。对于向量X，mean （X） 得到它的元素的算术平均值;
对于矩阵,mean （X）得到X各列元素的算术平均值，返回一个行向量。
nanmean 求忽略NaN的随机变量的算术平均值。 geomean 求随机变量的几何平均值。 harmmean 求随机变量的和谐平均值。 trimmean 求随机变量的调和平均值.
8.1.2 数据比较
在给定的一组数据中，还常要对它们进行最大、最小、中值的查找或对它们排序等操作。Mtalab中也有这样的功能函数。
max 求随机变量的最大值元素。
nanmax 求随机变量的忽略NaN的最大值元素。 min 求随机变量的最小值元素.
nanmin 求随机变量的忽略NaN的最小值元素。 median 求随机变量的中值.
nanmedian 求随机变量的忽略NaN的中值。 mad 求随机变量的绝对差分平均值。 sort 对随机变量由小到大排序. sortrows 对随机矩阵按首行进行排序。
range 求随机变量的值的范围,即最大值与最小值的差（极差）。
8.1.3 累和与累积
求向量或矩阵的元素累和或累积运算是比较常用的两类运算，在Matlab中可由以下函数实现。
sum 若X为向量，sum （X)为X中各元素之和，返回一个数值；若X为矩阵，
sum （X）为X中各列元素之和,返回一个行向量.
nansum 忽略NaN求向量或矩阵元素的累和。
cumsum 求当前元素与所有前面位置的元素和。返回与X同维的向量或矩阵。 cumtrapz 梯形累和函数。 prod 若X为向量，prod (X）为X中各元素之积，返回一个数值；若X为矩阵,prod
（X)为X中各列元素之积，返回一个行向量。
cumprod 求当前元素与所有前面位置的元素之积.返回与X同维的向量或矩阵。
8.2 离散型随机变量的概率及概率分布
0
8.2.1 几个常见分布
8.2.1.1 二项分布
设随机变量X的分布律为:
P ｛X = k } = Cnk p k （1—p) n- k k = 0， 1, 2， …, n
其中:0< p 〈1，n为独立重复试验的总次数，k为n次重复试验中事件A发生的次数，p为每次试验事件A发生的概率.则称X服从二项分布，记为X~B (n, p）。 8。2.1。2 Poisson分布
设随机变量X的分布律为：
P{xk}kk!ek0,1,2,;0
则称X服从参数为的Poisson分布，记为X~P （)。Poisson是二项分布的极限分布，当二项分布中的n较大而p又较小时,常用Poisson，= np。 8.2.1。3 超几何分布
设一批同类产品共N件，其中有M件次品，从中任取n （n≤N）件，其次品数X恰为k件的概率分布为：
knkCMCNM P{Xk},nCNk0,1,2,,min(n,M)
则称次品数X服从参数为 (N， M， n）的超几何分布。超几何分布用于无放回抽样，当N很大而n较小时,次品率p其中二项分布的p
M在抽取前后差异很小，就用二项分布近似代替超几何分布，NM.而且在一定条件下，也可用Poisson分布近似代替超几何分布。 N8。2.2 概率密度函数值
无论是离散分布还是连续分布，在Matlab中，都用通用函数pdf或专用函数来求概率密度函数值。而对于离散型随机变量，取值是有限个或可数个，因此，其概率密度函数值就是某个特定值的概率,即利用函数pdf求输入分布的概率。 1。 通用概率密度函数pdf计算特定值的概率
命令：pdf
格式为：Y = pdf (‘name'， k, A)
Y = pdf （‘name’， k, A, B) Y = pdf (‘name’, k， A, B, C）
说明：返回以name为分布,在随机变量X = k处，参数为A、B、C的概率密度值；对离散型随机变量X，返回X = k处的概率值，name为分布函数名。
常见的分布有：
name = bino(二项分布)，hyge(超几何分布），geo（几何分布），poiss（Poisson分布）。 2. 专用概率密度函数计算特定值的概率 （1）二项分布的概率值
命令：binopdf
格式:binopdf (k,n，p)
说明：等同于pdf (‘bino'， k, n, p）。n—试验总次数；p—每次试验事件A发生的概率；k-事件A发生k次。
（2）Poisson分布的概率值 命令:poisspdf
格式:poisspdf （k， Lambda)
说明：等同于pdf (‘poiss’, k, Lambda)，参数Lambda = np. (3）超几何分布的概率值 命令：hygepdf
1
格式:hygepdf （k， N, M, n)
说明:等同于pdf (‘hyge'， k， N， M， n)，N-产品总数,M—次品总数,n—抽取总数（n≤N)，k—抽得次品数.
3。 通用函数cdf用来计算随机变量X≤k的概率之和(累积概率值)
命令：cdf
格式：cdf (‘name’， k， A） cdf (‘name’, k, A, B） cdf (‘name’， k, A, B, C)
说明：返回以name为分布、随机变量X≤k的概率之和（即累积概率值），name为分布函数名。
4。 专用函数计算累积概率值（随机变量X≤k的概率之和，即分布函数)
（1）二项分布的累积概率值 命令：binocdf
格式：binocdf (k, n， p)
（2)Poisson分布的累积概率值 命令：poisscdf
格式：poisscdf （k, Lambda) （3）超几何分布的累积概率值 命令:hygecdf
格式：hygecdf （k, N, M， n)
8。2。2.1 二项分布
1。 求n次独立重复试验中事件A恰好发生k次的概率P。 命令：pdf 或 binopdf
格式：pdf (‘bino’， k， n, p) 或 binopdf （k， n, p）
说明:该命令的功能是计算二项分布中事件A恰好发生k次的概率。pdf为通用函数，bino表示二项分布，binopdf为专用函数，n为试验总次数，k为n次试验中，事件A发生的次数,p为每次试验事件A发生的概率。
2。 在n次独立重复试验中,事件A至少发生k次的概率P_s. 命令：cdf 或 binocdf
格式：cdf （‘bino’, k, n, p） 或 binocdf (k， n， p)
说明:该命令的功能是返回随机变量X≤k的概率之和（即累积概率值）。其中cdf为通用函数,binocdf为专用函数，n为试验总次数，k为n次试验中，事件A发生的次数，p为每次试验事件A发生的概率。
所以，至少发生k次的概率为
P_s = 1— cdf (‘bino’， k—1， n， p) 或 P_s = 1- binocdf (k—1, n， p） 例8-1 某机床出次品的概率为0。01，求生产100件产品中： （1）恰有1件次品的概率； (2）至少有1件次品的概率.
解:此问题可看作是100次独立重复试验,每次试验出次品的概率为0。01. (1）恰有1件次品的概率 在Matlab命令窗口键入：
>> p=pdf(’bino’,1，100，0.01） ％利用通用函数计算恰好发生k次的概率
p =
0。3697
2
或在Matlab命令窗口键入: 〉〉 p=binopdf（1，100，0.01) ％利用专用函数计算恰好发生k次的概率
p =
0。3697
(2)至少有1件次品的概率 在Matlab命令窗口键入：
〉> p=1-cdf（’bino’，0，100，0.01) % cdf是用来计算X≤k的累积概率值的通用函数，这里是计算X≥1的概率值。
p =
0。6340
或在Matlab命令窗口键入： 〉> p=1-binocdf（0,100,0。01) p =
0.6340
8。2.2。2 Poisson分布
在二项分布中，当n的值很大,p的值很小，而np又较适中时,用Poisson分布来近似二项分布较好（一般要求= np<10）。
1. n次独立重复试验中，事件A恰好发生k次的概率P_k。 命令:pdf
或 poisspdf
格式：pdf (‘poiss’, k, Lambda) 或 poisspdf （k, Lambda)
说明：在Matlab中，poiss表示Poisson分布。该命令返回事件恰好发生k次的概率. 2. n次独立重复试验中,事件A至少发生k次的概率P （1）累积概率值 命令：cdf
或 poisscdf
格式：cdf （‘poiss’， k， Lambda) 或 poisscdf （k， Lambda)
说明：该函数返回随机变量X≤k的概率之和，Lambda = np (2）A至少发生k次的概率P_k
P_k = 1- cdf （‘poiss’, k—1， Lambda) 或 P_k = 1- poisscdf （k-1， Lambda)
例8-2 自1875年到1955年中的某63年间，某城市夏季（5—9月间）共发生暴雨180次，试求在一个夏季中发生k次（k = 0, 1, 2, …， 8）暴雨的概率P k（设每次暴雨以1天计算）。
解：一年夏天共有天数为
n = 31+30+31+31+30 = 153
故可知夏天每天发生暴雨的概率约为p分布近似= np =
180，很小，n = 153较大,可用Poisson
63153180. 63在Matlab编辑器中编写M文件：LX0802.m p=input（’input p=’) n=input（’input n=') lambda=n*p
for k=1：9 ％循环变量的最小取值是从k = 1开始。 p_k（k)=poisspdf(k—1，lambda)； end
3
p_k
在Matlab的命令窗口键入LX0802，回车后按提示输入p和n的值,显示如下: input p=180/(63＊153） p =
0.0187 input n=153 n = 153 lamda =
2。8571 p_k =
Columns 1 through 7
0。0574 0.1641 0.2344 0。2233 0。1595 0。0911 0。0434 Columns 8 through 9 0。0177 0。0063
注意:在Matlab中，p_k （0)被认为非法，因此应避免。
例8—3 某市公安局在长度为t的时间间隔内收到的呼叫次数服从参数为t/2的Poisson分布，且与时间间隔的起点无关（时间以小时计）。
求：(1）在某一天中午12时至下午3时没有收到呼叫的概率； （2）某一天中午12时至下午5时至少收到1次呼叫的概率。 解：在此题中，Lamda = t/2
设呼叫次数X为随机变量，则该问题转化为： （1）求P｛X = 0｝； （2）求1—P{X≤0}。
解法一：在Matlab命令窗口键入：
〉> poisscdf (0，1.5) %X = 0表示0次呼叫，Lambda = t/2 = 1.5 ans =
0。2231
即（1)中没有收到呼叫的概率为0。2231. 〉〉 1—poisscdf (0，2。5) ans =
0.9179
即（2)中至少收到1次呼叫的概率为0。9179。 解法二:
由于呼叫次数X≤0就是呼叫0次，即X = 0。因此，此题也可用poisspdf求解。即： poisspdf （0, 1。5)和1—poisspdf (0， 2。5)。
8.2。2。3 超几何分布
1。 设N为产品总数，M为其中次品总数，n为随机抽取件数（n≤N），则次品数X恰为k件的概率p_k（k = 0， 1, 2, …, min (n， M））可由下列命令求得：
命令：pdf 或 hygepdf
格式：pdf (‘hyge’, k， N， M， n) 或 hygepdf （k， N, M, n） 2. 累积概率值的求法： 命令：cdf 或 hygecdf
格式：cdf （‘hyge', k， N， M, n） 或 hygecdf (k, N, M， n）
说明:该函数的功能是返回次品数X≤k的概率之和。
例8—4 设盒中有5件同样的产品,其中3件正品,2件次品,从中任取3件，求不能取得次品的概率。
4
解：在Matlab编辑器中编辑M文件:LX0802.m N=input('input N=’） M=input（’input M=’） n=input（'input n=’) for k=1:M+1
p_k=hygepdf（k—1,N,M，n) end
在Matlab的命令窗口键入LX0804，回车后按提示输入N、M、n的值，显示如下: input N=5 N = 5 input M=2 M = 2 input n=3 n = 3 p_k =
0.1000 p_k =
0。6000 p_k =
0。3000
这里，p_k=（0.1000 0。6000 0.3000）表示取到次品数分别为X = 0， 1， 2的概率.
8。3 连续型随机变量的概率及其分布
8.3.1 几个常见的分布
8。3。1。1 均匀分布
若随机变量X的概率密度为
1p(x)ba0axb其它，
则称X在区间［a, b]上服从均匀分布，记为X~U （a，
b）。
8。3。1。2 指数分布
若随机变量X的概率密度为
exp(x)0x0x0，其中0为常数
则称X服从参数为的指数分布,记为E （）。 8。3.1。3 正态分布
若随机变量X的概率密度为
22的正态分布，(0)是两个常数，其中，则称X服从参数为，记为X～N (,)。
8.3。1.4 分布
1xxex0 p(x)()， 其中0，0
0x0记为：(,)
5
p(x)1e(x)222x
8.3。1。5 分布
()1x(1x)1 p(x)()()0记为:(,)
8。3。1.6 2分布（卡方分布）
X的分布密度为
nx11x2e2n p(x)22(n)20n的2分布,记为：X～（n）。
8.3。1.7 t分布
若随机变量t的分布密度为
20x1其它， 其中0,0
x0， n为正整数，则称X为服从自由度为
x0
n1)1t2n22p(t)(1)nnn()2(t ，,
n为正整数，则称t为服从自由度为n的t分布,记为:T~ t (n）。 8。3.1.8 F分布
若随机变量x的分布密度为
n1n2n1n2n11)n12n22(2x2x0,n，n为正整数，则称X服从第n1n2p(x)12n1n2()()(n1xn2)222x00一自由度为n1,第二自由度为n2的F分布，记为：X~F （n1, n2)。
8。3。2 概率密度函数值
连续型随机变量：如果存在一非负可积函数p(x)0，使对于任意实数ab，X在区
P{aXb}p(x)dx,则函数p(x)称作随机变量X的概率密度间(a,b)上取值的概率为：
ab函数。通用函数pdf和专用函数用来求密度函数p(x)在某个点x处的值。
1。 利用概率密度函数值通用函数pdf计算 格式：pdf （‘name', x， A) pdf (‘name’， x, A, B）
pdf （‘name’, x， A, B, C）
说明：返回以name为分布的随机变量在X = x处、参数为A、B、C的概率密度函数值。name取值如下表8。1所示。
表8.1 常见通用函数密度函数表 name 分 布 unif 均匀分布密度函数 exp 指数分布密度函数 norm 正态分布密度函数 chi2 卡方（2）分布
6
logn 对数分布 nbin 负二项分布 ncf 非中心F分布 用函数计算概率nct 非中心t分布 ncx2 非中心卡方分布 表8。2 专用函数概率密度函数表 瑞利分布 函 数 名 调 用 形rayl 式 注 释 weib weibull( unifpdf unifpdf (x, a, b) ［a， b]韦伯）分布上均匀分布概率密度在X = x处的函数值 exppdf exppdf （x, Lambda) 指数分布概率密度在X = x处的函数值 normpdf normpdf (x, mu, sigma) 正态分布概率密度在X = x处的函数值 chi2pdf chi2pdf （x， n） 卡方分布概率密度在X = x处的函数值 tpdf tpdf (x, n） t分布概率密度在X = x处的函数值 fpdf fpdf (x, n1, n2) F分布概率密度在X = x处的函数值 gampdf gampdf (x, a, b) 分布概率密度在X = x处的函数值 betapdf betapdf （x， a， b) 分布概率密度在X = x处的函数值 lognpdf lognpdf (x, mu, sigma) 对数分布概率密度在X = x处的函数值 nbinpdf nbinpdf （x, R, P) 负二项分布概率密度在X = x处的函数值 ncfpdf ncfpdf （x, n1， n2, delta) 非中心F分布概率密度在X = x处的函数值 nctpdf nctpdf （x, n， delta) 非中心t分布概率密度在X = x处的函数值 ncx2pdf ncx2pdf （x, n, delta) 非中心卡方分布概率密度在X = x处的函数值 raylpdf raylpdf （x， b） 瑞利分布概率密度在X = x处的函数值 weibpdf weibpdf （x， a, b) weibull分布概率密度在X = x处的函数值 密度函数值，如表8.2所示。
例8—5 计算正态分布N （0,1）在点0.7733的概率密度值. 解：在Matlab命令窗口键入：
〉> pdf(’norm’,0。7733，0，1） ％利用通用函数 ans =
0.2958 〉〉 normpdf(0.7733，0，1) %利用专用函数 ans =
0。2958
两者计算结果完全相同。
例8-6 绘制卡方分布密度函数在n分别等于1，5，15时的图形。 解：在Matlab编辑器中编辑M文件：LX0806。m x=0：0。1:30； y1=chi2pdf(x，1); plot（x,y1，':'） hold on
y2=chi2pdf(x，5)； plot（x，y2,’+’） y3=chi2pdf（x,15）; plot（x，y3,'o’)
axis（[0，30,0,0.2］) xlabel（’图8。1’）
在命令窗口键入LX0806，回车后得结果为图8.1。
7
2. 利用专t或T f和F gam beta t分布 F分布密度函数 分布密度函数 分布
8.3.3 累积概率函数值（分布函数）
连续型随机变量的累积概率函数值是指随机变量X≤x的概率之和。即：
P{X≤x}=xp(t)dt
也就是连续型随机变量的分布函数F （x)，F (x）既可以用通用函数，也可用专用函数来计算.通常用这些函数计算随机变量落在某个区间上的概率和随机变量X的分布函数F （x)。 1。 利用通用函数cdf计算累积概率值
格式:cdf (‘name', x, A)
cdf (‘name’， x, A， B）
cdf （‘name’， x， A， B, C）
说明：返回随机变量X≤x的概率之和.name为上述分布函数名。 2。 利用专用函数计算累积概率值
其命令函数是在上述分布后面加上cdf，其用法同专用函数计算概率密度函数值。 如正态分布的累积概率值:
命令函数为：normcdf （x, mu， sigma）
则显示结果为 F （x） =xp(t)dt的值.
例8-7某公共汽车站从上午7：00起每15分钟来一班车。若某乘客在7：00到7:30间的任何时刻到达此站是等可能的，试求他侯车的时间不到5分钟的概率。
解：设乘客7点过X分钟到达此站,则X在[0,30］内服从均匀分布,当且仅当他在时间间隔（7：10，7:15）或（7:25，7：30）内到达车站时，侯车时间不到5分钟。故其概率为：
p = P{10 p1=unifcdf（15，0,30）—unifcdf(10,0，30）； p2=unifcdf(30,0，30）-unifcdf(25，0，30）； p=p1+p2 运行结果为： p = 1/3 例8-8 设X~N （3， 22),求P｛2〈X〈5｝，P{—4〈X〈10}，P{|X|>2｝，P{X〉3｝。 解：在Matlab编辑器中编辑M文件LX0808.m如下： 8 p1=normcdf(5，3,2）—normcdf（2，3,2） p2=normcdf（10,3，2)-normcdf（—4，3，2) p3=1—normcdf（2,3,2）+normcdf（-2，3,2) p4=1-normcdf（3,3,2） 运行结果为： p1 = 0。5328 p2 = 0.9995 p3 = 0。6977 p4 = 0.5000 例8-9 设随机变量X的概率密度为 c|x|1 p(x)1x2 0|x|1（1）确定常数c; 1122(3）求X的分布函数F (x)。 解：(1）在Matlab编辑器中建立M文件LX08091.m如下： syms c x p_x=c/sqrt（1—x^2)； F_x=int(p_x,x,—1,1) 运行结果为： F_x = pi*c 由pi＊c=1得 c=1/pi （2)在Matlab编辑器中建立M文件LX08092。m如下： syms x c='1/pi’; %’1/pi’不加单引号“'’”，其结果的表达式有变化. p_x=c/sqrt(1-x^2); format rat p1=int（p_x，x，—1/2，1/2) 运行结果为： p1 = 1/3 （3）在Matlab编辑器中建立M文件LX08093.m如下： syms x t c=’1/pi'; p_t=c/sqrt（1—t^2）； F_x=int（p_t，t，-1,x） 运行结果为: F_x = 1/2＊(2＊asin（x)+pi)/pi >> simple（F_x) ans = asin（x)/pi+1/2 （2）求X落在区间(,)内的概率； 9 x00arcsinx11x1 所以X的分布函数为：F （x） = 2x111例8-10 设lnX ～ N （1, 22），求P{X2}。 2解：利用对数分布累积专用函数，在Matlab命令窗口键入： >> p=logncdf(2，1，2）-logncdf（1/2，1，2） p = 0。2404 8.3.4 逆累积概率值 已知分布和分布中的一点，求此点处的概率值要用到累积概率函数cdf，当已知概率值而需要求对应概率的分布点时，就要用到逆累积概率函数icdf，icdf返回某给定概率值下随机变量X的临界值，实际上就是cdf的逆函数,在假设检验中经常用到。 即：已知F (x） = P{X≤x}，求x 逆累积概率值的计算有下面两种方法. 8。3。4。1 通用函数icdf 格式：icdf (‘name’， p, a1, a2， a3) 说明：返回分布为name，参数为a1, a2， a3累积概率值为p的临界值，这里name与前面相同。 如：p = cdf （‘name’, x， a1, a2, a3） 则:x = icdf (‘name’, p， a1， a2， a3) 例8-11设X~N （3, 22)，确定c,使得P{X>c} = P{X〈c}。 解：若要P{X>c} = P{X〈c｝,只需P｛X〉c｝ = P｛X〈c｝= 0。5 在Matlab命令窗口键入： 〉> c=icdf（'norm'，0.5，3,2) 运行结果为： c = 3 例8—12 在假设检验中，求临界值问题: 已知：0.05，查自由度为10的双边界检验t分布临界值。 解：在Matlab命令窗口键入： >> t0=icdf（'t'，0.025，10) 运行结果为: t0 = —2。2281 8.3。4。2 专用函数—inv 如：norminv （p， mu, sigma） ％正态逆累积分布函数，返回临界值.用法与前面类似。 关于常用临界值函数可查表8。3： 表8.3 常用临界值函数表 函 数 名 调 用 形 式 注 释 [a， b］上均匀分布逆累积分布函数，X为临unifinv x = unifinv (p, a, b) expinv norminv chi2inv tinv finv x = expinv (p, lambda） x = norminv （p， mu， sigma） x = chi2inv (p, n） x = tinv (p， n) x = finv （p， n1， n2) 界值 指数逆累积分布函数 正态逆累积分布函数 卡方逆累积分布函数 T分布逆累积分布函数 F分布逆累积分布函数 10 gaminv betainv logninv ncfinv 非中心t分布逆累积分布函数 nctinv 非中心卡方逆累积分布函数 ncx2inv 瑞利逆累积分布函数 raylinv 韦伯逆累积分布函数 weibinv 例8-13 公共汽车门的高度是按成年男子与车门顶碰头的机会不超过1%设计的。设男子身高X(单位：cm）～ N (175, 36)，求车门的最低高度。 解：设h为车门高度，X为男子身高，求满足条件P{X>h｝≤0。01的h，即P{X 188。9581 例8-14 设二维随机向量(X, Y)的联合密度为 x = gaminv （p， a, b) x = betainv (p， a, b) x = logninv (p, mu, sigma) x = ncfinv (p， n1， n2， delta) x = nctinv (p, n， delta) x = ncx2inv (p， n， delta) x = raylinv （p, b） x = weibinv （p, a, b） 分布逆累积分布函数 分布逆累积分布函数 对数逆累积分布函数 非中心F分布逆累积分布函数 e(xy)x0,y0 p(x,y) 其它0求:（1）P{0〈X〈1， 0 f=exp（-x-y); p_XY=int（int(f,y,0，1)，x,0,1) p_G=int（int(f,y，0，1—x)，x,0，1) 运行结果为： p_XY = exp(—2)-2*exp(-1)+1 p_G = —2*exp(-1)+1 8。4 数字特征 随机变量的数字特征是概率统计学的重要内容。在对随机变量的研究中，如果对随机变量的分布不需要作全面的了解，那么只需要知道它在某一方面的特征就够了.这些特征的数就是随机变量的数字特征。 8.4。1 随机变量的期望 期望是随机变量的所有可能取值乘以相应的概率值之和，即 EXxkpk k1其中pk是对应于xk的概率,即权重。 1是常用的情况：给定一组样本值x = ［x1, x2， …, xn］ n1 EXxk nk1此时，期望称为样本均值. 8。4.1.1 离散型随机变量X的期望计算 pk 11 1. 函数：sum ％求和 格式：sum （X) 说明：若X为向量，则sum （X）为X中的各元素之和,返回一个数值； 若X为矩阵，则sum （X)为X中各列元素之和，返回一个行向量。 2. 函数：mean 格式:mean （X) 说明:若X为向量,则mean (X）为X中的各元素的算术平均值，返回一个数值； 若X为矩阵,则mean （X)为X中各列元素的算术平均值，返回一个行向量。 例8-15设随机变量X的分布律为： X -2 -1 0 1 p 0.3 0。1 0。2 0。1 求EX，E (X2-1）。 解：在Matlab编辑器中建立M文件LX0815。m: X=［-2 —1 0 1 2]； p=［0.3 0.1 0.2 0.1 0。3］; EX=sum（X＊p'） Y=X。^2-1; EY=sum(Y＊p') 运行结果为： EX = 0 EY = 1.6000 或 X=［—2 —1 0 1 2]； p=[0。3 0.1 0.2 0。1 0。3］； EX=sum(X。*p) Y=X。^2—1； EY=sum（Y.*p） EX = 0 EY = 1。6000 例8-16 随机抽取6个滚珠测得直径如下:（单位：mm) 14。70 15.21 14。90 14.91 15。32 15。32 试求样本均值。 解：在Matlab命令窗口键入： >〉 X=[14。70 15。21 14。90 14.91 15.32 15.32]； >〉 mean（X） ans = 15.0600 8。4.1。2 连续型随机变量的期望 若随机变量X的概率密度为p (x)，则X的期望为: EXxp(x)dx 若下式右端积分绝对收敛，则X的函数f (X）的期望为： Ef(X)f(x)p(x)dx 例8—17 已知随机变量X的概率密度 p(x)3x2 0x10其它 12 2 0。3 求EX和E (4X-1）. 解：在Matlab编辑器中建立M文件LX0817。m： syms x p_x=3*x^2； EX=int(x＊p_x，0,1) EY=int((4*x—1)*p_x，0，1) 运行结果为： EX = 3/4 EY = 2 例8-18 设随机变量X的概率密度为 p(x)1|x|e,2x 求EX. 解：在Matlab编辑器中建立M文件LX0818。m： syms x p_x=1/2*exp(—abs（x))； EX=int（x*p_x，—inf，inf） 运行结果为： EX = 0 8.4。2 方差 方差是随机变量的个别偏差的平方的期望： DX = E (X-EX）2 = EX2—(EX）2 标准差:(X)DX 对于样本x = [x1， x2, …， xn］，有 1n2(xix)2 样本方差：Sn1i11n样本标准差：S(xix)2 n1i18.4.2.1 离散型随机变量的方差 1。 方差 DX = E (X—EX）2 = EX2-（EX）2 在Matlab中用sum函数计算： 设X的分布律为P{X = xk} = pk,k = 1， 2, … 则方差 DX = sum （X - EX)。^2.＊p 或 DX = sum （X 。^2。＊p) - (EX).^2 标准差：(X)DX = sqrt (DX） 例8-19 设随机变量X的分布律为 X -2 -1 0 p 0.3 0.2 0。1 求DX，D (X2-1)。 解：在Matlab编辑器中建立M文件LX0819.m： X=[-2 -1 0 1 2]; p=［0.3 0。1 0.2 0。1 0。3]; EX=sum(X.*p） Y=X.^2—1 EY=sum（Y。*p) DX=sum（X.^2.＊p)—EX。^2 13 1 0.1 2 0.3 DY=sum（Y.^2.*p）-EY.^2 运行结果为： EX = 0 Y = 3 0 -1 0 3 EY = 1。6000 DX = 2。6000 DY = 3。0400 2. 样本方差 设随机变量X的样本为x = [x1, x2, …, xn],由于X取xi的概率相同且均为1/n，因此可以用上面的方法计算方差。另一方面，在Matlab中又有专门的函数var计算样本方差。 函数:var ％计算一组采集数据即样本的方差 1n2(xix)2，若X为向量，则返回向量格式：var (X) ％var (X) = Sn1i1的样本方差；若X为矩阵,则返回矩阵列向量的样本方差构成的行向量。 var （X, 1） ％返回向量（矩阵)X的简单方差（即置前因子为1/n的方差） var （X, w) ％返回向量（矩阵）X的以w为权重的方差 函数：std %计算一组采集数据即样本的标准差 格式:std (X） ％返回向量（矩阵）X的样本标准差即: 1nstd （X) = S(xix)2 n1i1std (X, 1） ％返回向量(矩阵）X的标准差（置前因子为1/n) std （X, 0) %与std (X）相同 std （X， flag, dim) %返回向量（矩阵)X中维数为dim的标准差值，其中flag=0 时,置前因子为1/（n—1);否则置前因子为1/n. 例8-20 求下列样本的样本方差和样本标准差，方差和标准差 14.70 15.21 14.90 14.91 15。32 15.32 解：在Matlab编辑器中建立M文件LX0820。m： X=［14.70 15。21 14。90 14.90 15.32 15。32］; DX=var（X,1) ％方差 sigma=std（X，1) %标准差 DX1=var（X） %样本方差 sigma1=std（X） ％样本标准差 运行结果为： DX = 0.0559 sigma = 0。2364 DX1 = 0。0671 sigma1 = 0.2590 8.4.2。2 连续型随机变量的方差 利用 DX = E (X—EX)2 = EX2-（EX)2 求解。 设X的概率密度为p （x) 14 则 EX DX或 DXxp(x)dx (xEX)2p(x)dx xp(x)dx(2xp(x)dx)2 在Matlab中，视具体情况实现. 例8—21 设X的密度函数为 1 p(x)1x20|x|1|x|1 求DX，D （2X+1)。 解:在Matlab编辑器中建立M文件LX0821.m： syms x px=1/（pi*sqrt(1-x^2））； EX=int（x＊px,-1，1) DX=int(x^2*px,-1，1)-EX^2 y=2＊x+1; EY=int(y＊px，—1，1) DY=int（y^2*px,—1，1）—EY^2 运行结果为: EX = 0 DX = 1/2 EY = 1 DY = 2 8.4.3 常用分布的期望与方差求法 在统计工具箱中，用‘stat’结尾的函数可以计算给定参数的某种分布的均值和方差.见下表： 表8。4 期望和方差表 函数名 调 用 形 式 参 数 说 明 函 数 注 释 betastat binostat chi2stat expstat fstat gamstat geostat hygestat lognstat poisstat normstat tstat unifstat [M, V]= betastat (A, B） [M， V]=binostat (N， p) M为期望值,V为方差值； A,B为分布参数 分布的期望与方差 N为试验次数； p为二项分布概率 [M， V]= chi2stat (nu) nu为卡方分布参数 ［M， V]= expstat (mu） mu为指数分布参数 [M, V]= fstat （n1， n2) n1, n2为F分布的两个自由度 [M， V]= gamstat （A， B) A， B为分布的参数 [M, V]= geostat （p） p为几何分布的几何概率参数 ［M, V]= hygestat (M，K,N) M, K, N为超几何分布的参数 mu为对数分布的均值； [M， V]= lognstat (mu, sigma) sigma为标准差 ［M, V]=poisstat (lambda) lambda为Poisson分布的参数 mu为正态分布的均值; [M, V]= normstat （mu, sigma） sigma为标准差 nu为t分布的参数 [M, V]= tstat (nu) a, b为均匀分布的分布区间端［M, V]= unifstat (a， b） 点值 二项分布的期望与方差 卡方分布的期望与方差 指数分布的期望与方差 F分布的期望与方差 分布的期望与方差 几何分布的期望与方差 超几何分布的期望与方差 对数分布的期望与方差 Poisson分布的期望与方差 正态分布的期望与方差 t分布的期望与方差 均匀分布的期望与方差 有了上面的表格，各函数的用法也就一目了然了。下面举几个例子. 15 例8-22 求参数为0。12和0。34的分布的期望和方差。 解：在Matlab命令窗口键入： >> ［m,v]=betastat（0。12，0.34) 结果为： m = 0.2609 v = 0。1321 例8—23 按规定,某型号的电子元件的使用寿命超过1500小时为一级品，已知一样品20只，一级品率为0。2。问这样品中一级品元件的期望和方差为多少? 解：分析可知此电子元件中一级品元件分布为二项分布,可使用binostat函数求解。在Matlab命令窗口键入： 〉〉 [m,v］=binostat（20,.2） 结果为： m = 4 v = 3.2000 结果说明一级品元件的期望为4，方差为3.2. 例8—24 求参数为8的Poisson分布的期望和方差。 解： >> [m，v]=poisstat(8) m = 8 v = 8 由此可见Poisson分布参数的值与它的期望和方差是相同的。 8。5 二维随机向量的数字特征 8。5.1 期望 1. 若(X, Y)的联合分布律为 P{X = xi, Y = yj ｝= pij i = 1, 2, …； j = 1， 2， … 则Z = f (X， Y）的期望为： EZ = E f （X， Y） = f(xi,yj)pij ij2. 若(X， Y）的联合密度为p （x, y）,则Z = f (X， Y)的期望为： EZ = E f （X, Y） = EX EYf(x,y)p(x,y)dxdy 3. 若（X， Y）的边缘概率密度为pX （x），pY （y)，则 xpX(x)dxxp(x,y)dxdy ypY(y)dyyp(x,y)dxdy (xEX) DX DY(xEX)2pX(x)dx2p(x,y)dxdy p(x,y)dxdy (yEY)2pY(y)dy(yEY)2例8—25 设（X， Y)的联合分布为 Y -1 X 1 2 16 -1 2 5 203 202 203 206 201 20Z = X-Y，求EZ。 解：在Matlab编辑器中建立M文件LX0825.m： X=［-1 2］; Y=[-1 1 2]； for i=1:2 for j=1：3 Z(i，j)=X（i）-Y（j)； end end ％该循环计算X—Y的值Z p=［5/20 2/20 6/20;3/20 3/20 1/20］； EZ=sum（sum（Z.＊p)) %将Z与p对应相乘相加 运行结果为： EZ = —0。5000 例8-26 射击试验中，在靶平面建立以靶心为原点的直角坐标系,设X、Y分别为弹着点的横坐标和纵坐标，它们相互独立且均服从N （0， 1)，求弹着点到靶心距离的均值。 解：弹着点到靶心的距离为Z其联合分布密度为 X2Y2，求EZ 112(x2y2) p(x,y)ex,y 2在Matlab编辑器中建立M文件LX0826。m： syms x y r t pxy=1/（2＊pi）*exp（-1/2*（x。^2+y。^2））； EZ=int（int(r＊1/(2*pi）*exp(—1/2＊r^2)＊r,r，0，inf),t，0，2＊pi) %利用极坐标计算较简单 运行结果为： EZ = 1/2*2^（1/2）＊pi^（1/2） 2即 EZ = 28。5.2 协方差 对于二维随机向量（X, Y)，期望EX、EY分别反映X、Y各自的均值，而方差DX、DY也仅仅反映分量X、Y对各自均值的离散程度.因此还需要研究X与Y之间相互联系的程度。协方差是体现这一程度的一个很重要的概念。 设（X， Y)是一个二维随机向量，若E[(X-EX）(Y—EY）]存在，则称之为X,Y的协方差，记为cov (X, Y）或XY。 即：cov （X， Y） = E[（X-EX）（Y—EY)] = E (XY)—EX*EY 特别地：cov (X， X) = E[(X-EX）2] = EX2 - (EX）2 cov （Y, Y) = E［（Y-EY)2] = EY2 — (EY)2 Matlab提供了求样本协方差的函数： cov (X) %X为向量时,返回此向量的方差;X为矩阵时，返回此矩阵的协方差矩阵， 此协方差矩阵对角线元素为X矩阵的列向量的方差值. cov （X, Y） %返回X与Y的协方差，且X与Y同维. cov （X, 0) ％返回X的样本协方差，置前因子为1/(n-1）与cov (X）相同。 cov （X， 1) %返回X的协方差,置前因子为1/n。 cov （X， Y）与cov (X, Y， 1)的区别同上。 17 说明：用命令函数cov时，X，Y分别为样本点。 例8—27 设（X， Y)的联合密度为 1(xy)0x2,0y2 p(x,y)8 其它0求DX，DY和cov (X， Y）。 解：EX EYxpX(x)dxypY(y)dyxp(x,y)dxdy yp(x,y)dxdy 在Matlab编辑器中建立M文件LX0827。m: syms x y pxy=1/8*(x+y)； EX=int（int(x＊pxy，y，0，2)，0，2） EY=int（int（y＊pxy，x,0,2),0，2） EXX=int(int（x^2*pxy,y,0，2）,0,2) EYY=int（int（y^2*pxy,x，0,2)，0,2) EXY=int(int(x*y＊pxy，x,0，2）,0,2） DX=EXX-EX^2 DY=EYY-EY^2 DXY=EXY—EX*EY 运行结果为： EX = 7/6 EY = 7/6 EXX = 5/3 EYY = 5/3 EXY = 4/3 DX = 11/36 DY = 11/36 DXY = -1/36 例8-28 求向量a = ［1 2 1 2 2 1］的协方差。 解：在Matlab命令窗口键入： 〉〉 a=［1 2 1 2 2 1］; 〉〉 cov(a) ans = 0.3000 例8-29 求矩阵的协方差. 解:在Matlab窗口键入： 〉> d=rand(2，6） 18 d = 0。9218 0。1763 0.9355 0。4103 0。0579 0。8132 0.7382 0.4057 0.9169 0。8936 0.3529 0.0099 〉> cov1=cov（d） cov1 = 0。0169 -0.0211 0.0017 -0.0444 —0。0271 0.0737 -0。0211 0.0263 —0。0021 0.0555 0。0338 —0。0922 0。0017 -0。0021 0.0002 -0.0045 —0。0027 0.0075 -0.0444 0.0555 —0。0045 0。1168 0。0713 —0。1942 —0。0271 0.0338 —0.0027 0。0713 0.0435 —0.1185 0。0737 -0.0922 0。0075 —0。1942 -0.1185 0。3226 8。5。3 相关系数 相关系数是体现随机变量X和Y相互联系程度的度量。 设（X， Y)的协方差为cov (X， Y)，且DX〉0，DY〉0，则称为X与Y的相关系数，记为XY. 当XY= 0时,称X与Y不相关. Matlab提供了求样本相关系数的函数。 corrcoef (X,Y） %返回列向量X,Y的相关系数。 corrceof (X） ％返回矩阵X的列向量的相关系数矩阵。 例8-30 设（X， Y)的联合分布律为 Y -1 1 2 X XYXXYY即 cov(X,Y)DXDY526 —1 202020331 2 202020求X与Y的协方差XY及相关系数XY。 解:在Matlab编辑器中建立M文件LX0830.m： format rat ％有理格式输出 X=[—1 2]; Y=[-1 1 2］； PXY=［5/20 2/20 6/20；3/20 3/20 1/20］; ％X、Y的联合分布 PX=sum(PXY’) %X的边缘分布 PY=sum(PXY) ％Y的边缘分布 EX=sum（X.＊PX) %X的期望 EY=sum（Y.＊PY) EXX=sum(X。^2。＊PX） ％计算EX2 EYY=sum(Y。^2.*PY） DX=EXX-EX^2 ％X的方差 DY=EYY—EY^2 XY=［1 -1 -2；—2 2 4]； ％XY的取值 EXY=sum(sum(XY.＊PXY）) DXY=EXY—EX*EY %X与Y的协方差 ro_XY=DXY/sqrt（DX＊DY） %X与Y的相关系数 运行结果为: PX = 19 13/20 7/20 PY = 2/5 1/4 7/20 EX = 1/20 EY = 11/20 EXX = 41/20 EYY = 41/20 DX = 819/400 DY = 699/400 EXY = -1/4 DXY = -111/400 ro_XY = -365/2488 即：XY= —111/400，XY= -365/2488 例8—31 设（X， Y)在单位圆G={（x，y)|x2+y2≤1}上服从均匀分布,即有联合密度 1x2y21 p(x,y) x2y210求XX,XY,YY,XY. 解:EXxp(x,y)dxdyrcosp(x,y)rdrd GG在Matlab编辑器中建立M文件LX0831。m： syms x y r t pxy='1/pi’； ％若’1/pi’不加单引号，其结果表达式将较繁 EX=int（int（r^2*cos（t)＊pxy,r，0，1）,0，2*pi） EY=int(int（r^2＊sin（t)*pxy，r,0,1）,0，2*pi) EXX=int（int(r^3*cos(t）^2*pxy，r,0，1），0,2＊pi） EYY=int（int（r^3*sin（t）^2*pxy，r,0,1）,0,2＊pi) EXY=int(int（r^3＊cos(t）＊sin(t）*pxy，r,0,1）,0,2*pi） DX=EXX—EX^2 DY=EYY—EY^2 DXY=EXY-EX*EY ro_XY=DXY/sqrt(DX*DY） 运行结果为： EX = 0 EY = 0 EXX = 1/4 EYY = 1/4 EXY = 0 DX = 20 1/4 DY = 1/4 DXY = 0 ro_XY = 0 8.6 统计直方图 函数 hist ％直角坐标系下的统计直方图 格式：hist （X， n) 说明：X为统计数据，n表示直方图的区间数，缺省值n =10。 函数 rose ％极坐标系下角度直方图 格式:rose (theta， n） 说明：n是在［0, 2]范围内所分区域数，缺省值n =20；theta为指定的弧度数据。 例8—32 某食品厂为加强质量管理，对生产的罐头重量X进行测试，在某天生产的罐头中抽取了100个，其重量测试数据记录如下： 342 340 348 346 343 342 346 341 344 348 346 346 340 344 342 344 345 340 344 344 343 344 342 343 345 339 350 337 345 349 336 348 344 345 332 342 342 340 350 343 347 340 344 353 340 340 356 346 345 346 340 339 342 352 342 350 348 344 350 335 340 338 345 345 349 336 342 338 343 343 341 347 341 347 344 339 347 348 343 347 346 344 345 350 341 338 343 339 343 346 342 339 343 350 341 346 341 345 344 342 试根据以上数据作出X的频率直方图。 解：在Matlab编辑器中建立M文件LX0832.m： X=［342 340 348 346 343 342 346 341 344 348 。.. 346 346 340 344 342 344 345 340 344 344 ..。 343 344 342 343 345 339 350 337 345 349 ..。 336 348 344 345 332 342 342 340 350 343 .。。 347 340 344 353 340 340 356 346 345 346 .。. 340 339 342 352 342 350 348 344 350 335 ..。 340 338 345 345 349 336 342 338 343 343 。。. 341 347 341 347 344 339 347 348 343 347 。.. 346 344 345 350 341 338 343 339 343 346 .。。 342 339 343 350 341 346 341 345 344 342］; hist（X,13） 结果为： 21 8.7 参数估计 通常，一个随机变量的分布可由某些参数决定，但在实际问题中要想知道一个分布的参数的精确值是很困难的，因此，需要对这些参数的取值作出估计.参数估计就是通过随机样本去估计参数的取值以及参数的取值范围. Matlab的统计工具箱中采用极大似然法给出了常用概率分布的参数的点估计和区间估计值。另外还提供了部分分布的对数似然函数的计算功能。 8。7.1 点估计 点估计是用样本去估计总体参数取值的大小，这里有两种常用方法：样本数字特征法和矩估计法。 8。7.1。1 样本数字特征法 1n用样本均值xxi作为总体均值EX的估计值; ni11n(xix)2作为总体方差DX的估计值。 用样本方差Sn1i1在Matlab中,样本x = [x1， x2,…， xn］，则 样本均值：mx = 1/n＊sum (x） 样本方差：S2 = 1/(n-1)＊sum (（x-mx).^2） 8。7.1.2 矩估计法 这里主要介绍总体为正态分布的参数,2的矩估计。 2设X~N (,2），X1, X2,…, Xn为其样本，则,2的矩估计量为： 1nˆxix ni11n (xix)2 ni1在Matlab中，样本x = ［x1, x2，…, xn］，则 样本均值：mx = 1/n*sum (x） 样本方差：sigma = 1/n＊sum ((x-mx).^2） 28.7.2 最大似然估计法 Matlab统计工具箱中给出了最大似然法估计常用概率分布的参数的点估计和区间估计值函数，还提供了部分分布的对数似然函数的计算功能。 8.7。2。1 常用分布的参数估计函数(见下表） 表8.5 参数似然估计函数表 函数名 调 用 形 式 函 数 说 明 binofit poissfit normfit betafit unifit expfit gamfit binofit （X, N) ［PHAT, PCI] = binofit （X， N， ALPHA) poissfit (X) [LAMBDAHAT, LAMBDACI］= poissfit (X，) normfit （X, ALPHA) ［MUHAT, SIGMAHAT, MUCI， SIGMACI] = normfit (X， ALPHA) betafit (X) [PHAT， PCI］ = betafit (X, ALPHA) unifit (X, ALPHA) [AHAT, BHAT, ACI, BCI] = unifit （X, ALPHA) expfit (X） [MUHAT， MUCI］ = expfit （X, ALPHA) gamfit （X） [PHAT, PCI］ = gamfit （X， ALPHA） 22 二项分布的最大似然估计 返回水平的参数估计和置信区间 泊松分布的最大似然估计 返回水平的参数和置信区间 正态分布的最大似然估计 返回水平的期望、方差和置信区间 分布的最大似然估计 返回最大似然估计值和水平的置信区间 均匀分布的最大似然估计 返回水平的参数估计和置信区间 指数分布的最大似然估计 返回水平的参数估计和置信区间 分布的最大似然估计 返回最大似然估计值和水平的置信区间 weibfit mle weibfit （DATA， ALPHA） [PHAT， PCI］ = weibfit (DATA， ALPHA） PHAT = mle (DIST, DATA) [PHAT, PCI］ = mle (DIST， DATA， ALPHA, PI） 韦伯分布的最大似然估计 返回水平的参数及其区间估计 DIST分布的最大似然估计 返回最大似然估计值和水平的置信区间 说明:·各函数返回已给数据向量的参数最大似然估计值和置信度为（1-）×100%的置信区间。的默认值为0。05，即置信度为95%； ·在mle函数中,参数dist可为各种分布函数名，可实现各分布的最大似然估计. 例如分布： 函数:betafit 功能：分布数据的参数a和b的最大似然估计值及其置信区间. 格式：PHAT = betafit （X) [PHAT， PCI］= betafit (X, ALPHA） 说明：PHAT为样本X的分布参数a和b估计值；PCI为样本X的分布参数a和b的置信区间，是一个2×2矩阵，其第1列为参数a的置信下界和上界，第2列为为参数b的置信下界和上界，ALPHA为显著水平，(1—)×100％为置信度。 函数：mle 功能：求分布参数的最大似然估计值 格式:phat = mle （‘dist’, X） ［phat, pci］ = mle (‘dist’, X） ［phat, pci］ = mle （‘dist'， X， alpha） ［phat, pci］ = mle （‘dist', X, alpha, pl） ％仅用于二项分布,pl为试验次数。 说明：dist可为各种分布函数名,如:beta (分布)、bino （二项分布),X为数据样本，alpha为显著水平，（1-)×100%为置信度. 例8—33 随机产生100个分布数据，相应的分布参数真值为4和3。求4和3的最大似然估计值和置信度为99%的置信区间。 解：在Matlab编辑器中建立M文件LX0833.m： X=betarnd(4，3,100，1) %随机产生100个分布数据,参数为4和3。 [phat，pci]=betafit（X，0.01) 运行结果为： X = %这些数据只有一列，这里为了节约版面而改为4列数据(使用矩 阵重置命令reshape (X， 25, 4）） 0.3658 0.6519 0.6081 0.9078 0。3643 0。8469 0.5968 0。4355 0.3699 0。5699 0。2124 0。4981 0.4072 0.2345 0。4700 0.5750 0。6875 0。6418 0。6228 0.5098 0.7258 0。7138 0.8713 0。7770 0。6845 0.3643 0。4154 0。7091 0.5608 0。5030 0。5983 0.5150 0.8293 0。6394 0.6324 0.4216 0.2735 0.3465 0。5696 0.5543 0。6139 0.5409 0.5737 0.5949 0.5499 0。6392 0.7139 0.4601 0。4019 0.6719 0.5702 0.4127 0.5287 0。5353 0。8848 0。5694 0。2029 0.5285 0.6796 0.5562 0.5193 0。7248 0。6908 0。7405 0。7569 0.8543 0。2363 0。6161 0.7796 0。4654 0。3605 0。7372 23 0.5012 0.6840 0。4441 0。1429 0.7392 0。6577 0.6327 0。3682 0.7025 0.4687 0.6471 0。7881 0.4492 0.7995 0。2292 0。7464 0。6360 0.5585 0.5740 0.5893 0.6985 0.4931 0。4393 0.5544 0。4263 0.6238 0.4507 0.6960 phat = 4。6613 3.5719 pci = 3。1123 2.3336 6.2103 4.8102 说明：数据4。6613和3。5719为参数4和3的估计值；pci的第1列为参数4的置信区间，第2列为参数3的置信区间。随机产生的数据不同,其估计值和置信区间就不一样。 例8-34 设某种油漆的9个样品，其干燥时间（以小时计）分别为 6。0 5。7 5。8 6.5 7。0 6。3 5.6 6。1 5。0 设干燥时间总体服从正态分布N （,2),求和的置信度为0。95的置信区间（未知)。 解：在Matlab命令窗口键入： >> X=[6.0 5。7 5.8 6。5 7。0 6。3 5。6 6.1 5。0]; 〉> [muhat，sigmahat，muci,sigmaci］=normfit(X，0。05） muhat = 6 %的最大似然估计值 sigmahat = 0。5745 %的最大似然估计值 muci = %的置信区间 5。5584 6。4416 sigmaci = %的置信区间 0。3880 1。1005 此解说明的最大似然估计值为6,置信区间为[5.5584，6。4416］;的最大似然估计值为0。5745,置信区间为[0。3880，1.1005］。 例8—35 分别使用金球和铂球测定引力常数 (1)用金球测定观察值为：6.683 6。681 6.676 6.678 6.679 6.672 (2）用铂球测定观察值为:6.661 6.661 6。667 6.667 6.664 设测定值总体服从N (,2），和为未知.对（1）(2)两种情况分别求和为的置信度为0.9的置信区间。 解：在Matlab命令窗口键入： 〉> j=［6。683 6.681 6.676 6.678 6。679 6.672]; 〉> b=[6。661 6.661 6.667 6。667 6。664]; 〉〉 ［muhat，sigmahat，muci，sigmaci］=normfit（j，0.1） ％金球测定的估计 muhat = 6.6782 sigmahat = 0。0039 muci = 6。6750 6。6813 sigmaci = 0。0026 0。0081 24 说明金球测定数据的置信度为0.9的和置信区间为 ：[6。6750，6.6813］ ：［0.0026,0。0081］ >> [muhat，sigmahat，muci,sigmaci］=normfit（b，0。1） %铂球测定的估计 muhat = 6。6640 sigmahat = 0.0030 muci = 6。6611 6.6669 sigmaci = 0.0019 0.0071 说明铂球测定数据的置信度为0.9的和置信区间为 ［6.6611,6.6669］ ： ：[0。0019，0.0071］ 8。7。2.2 对数似然函数 Matlab统计工具箱提供了分布,分布，正态分布和韦伯分布的负对数似然函数值的求取函数. 1. betalike 功能：分布负对数似然函数。 格式:logL = betalike (params, data) [logL， info] = betalike （params， data) 说明：logL = betalike （params, data)返回分布对数似然函数负值。其中params为包含分布的参数a， b的矢量［a， b］，data为服从分布的样本数据，logL的长度与数据data的长度相同. ［logL, info］ = betalike （params, data)则同时给出了Fisher信息矩阵info。Info的对角线元素为相应参数的渐进方差。 betalike是分布最大似然估计的实用函数.似然函数假设数据样本中，所有的元素相互独立。因为betalike返回负对数似然函数，用fminsearch函数最小化betalike与最大似然估计的功能是相同的。 示例： 〉>r = betarnd（4,3，100,1）； ％随机产生的分布数据 〉〉［logl，avar] = betalike（[3.9010 2.6193］，r） logl = -33。0514 avar = 0.2856 0。1528 0。1528 0。1142 2。 gemlike 功能：分布的负对数似然函数。 格式：logL = gamlike （params， data) ［logL, info］ = gamlike （params, data） 说明:logL = gamlike （params， data）返回由给定样本数据data确定的分布的参数params(即[a， b])的负对数似然函数值。logL的长度与数据data的长度相同. ［logL, info］ = gamlike (params, data）则同时给出了Fisher信息矩阵info。Info的对角线元素为相应参数的渐进方差。 gamlike是分布的最大似然估计工具函数。因为gamlike返回负对数似然函数值,故用用fminsearch函数将gamlike最小后，其结果与最大似然估计是相同的。 25 示例: >>a = 2； b = 3; >>r = gamrnd(a，b，100，1）； 〉〉［logL，info］ = gamlike（[2.1990 2。8069]，r) logL = 267.5585 info = 0.0690 —0。0790 -0.0790 0。1220 3. normlike 功能:正态分布的负对数似然函数。 格式：logL = normlike (params, data） 说明：与betalike和gamlike的功能类似，不再赘述。params参数中，params （1）为正态分布的参数mu，params （2)为参数sigma。 4. weiblike 功能：weibull分布的负对数似然函数。 格式：logL = weiblike （params， data) [logL, info] = weiblike (params， data) 说明：weibull分布的负对数似然函数定义为： logLlogf(a,b|xi)logf(a,b|xi) i1i1nnlogL = weiblike （params, data）返回参数为a, b的weibull分布在给定数据点data时的负对数似然函数值。其中params (1) = a，params （2) = b。 ［logL， info] = weiblike (params, data)则增加了Fisher信息矩阵info。info的对角线元素为相应参数的渐进方差. 示例: >>r = weibrnd（0。5，0.8,100,1); 〉〉[logL，info] = weiblike（［0.4746 0。7832],r) logL = 203。8216 info = 0。0021 0.0022 0.0022 0。0056 8.8 假设检验 在Matlab中,假设检验问题都提出两种假设：即原假设和备择假设.对于正态总体均值的假设检验给出了检验函数: ztest 已知2，检验正态总体均值; ttest 未知2，检验正态总体均值； ttest2 两个正态总体均值比较. 对于一般连续型总体一致性的检验，给出了检验方法——秩和检验,由函数ranksum实现。 8。8。1 单个正态总体N （,2）的假设检验 8。8.1.1 2已知，对期望的假设检验——U检验法 函数:ztest 格式：H = ztest （X, m， sigma） H = ztest (X, m, sigma, alpha） ［H, sig， ci] = ztest （X, m, sigma， alpha， tail） 说明：X——样本 26 m—-期望值0 sigma—-正态总体标准差 alpha——检验水平（默认为0.05） tail--备选假设的选项，有三种情况： tail = 0（缺省）——m tail = 1——m tail = -1——m 即：tail = 0为双边检验，其余为单边检验问题。 H——检验结果，两种情况： H = 0-—在水平下,接受原假设，或假设相容; H = 1-—在水平下,拒绝原假设，或假设不相容。 sig——当原假设为真时（即m成立）,得到观察值的概率,当sig为小概率时， 则对原假设提出质疑 ci——均值的置信度为1的置信区间 例8-36 某车间用一台包装机包装葡萄糖，包得的袋装糖重是一个随机变量，它服从正态分布.当机器正常时，其均值为0。5kg，标准差为0。015。某日开工后检验包装机是否正常,随机地抽取所包装的糖9袋，称得净重为：（kg) 0.497 0.506 0。518 0。524 0。498 0.511 0。52 0.515 0。512 问机器工作是否正常？ 解:（1）分析 总体和已知，则可设样本的= 0。015，于是X~N (， 0.0152),问题就化为根据样本值来判断=0。5还是≠0.5。为此，提出假设 原假设: H0:00.5 备择假设: H1：0 （2)Matlab实现 >> X=［0.497 0.506 0.518 0.524 0.498 0.511 0。52 0.515 0.512］； ％注意：此处数据X 只能为向量而非矩阵 〉〉 ［H,sig]=ztest(X，0。5，0。015,0。05，0） H = 1 sig = 0。0248 结果H =1，说明在0。05的水平下，可拒绝原假设，即认为这天包装机工作不正常。 8.8。1。2 2未知，对期望的假设检验—— t检验法 函数:ttest 格式：H = ttest （X， m, alpha） ［H, sig, ci］ = ttest (X， m, alpha, tail) 说明：X——样本 m—-期望值0 alpha——检验水平（默认为0。05） tail——备选假设的选项,有三种情况: tail = 0（缺省）—-m tail = 1-—m tail = -1--m 即：tail = 0为双边检验，其余为单边检验问题. H—-检验结果，两种情况： H = 0——在水平下，接受原假设，或假设相容； H = 1——在水平下,拒绝原假设，或假设不相容. 27 sig—-当原假设为真时（即m成立），得到观察值的概率,当sig为小概率时， 则对原假设提出质疑 ci--均值的置信度为1的置信区间 例8—37 某种电子元件的寿命X（以小时计）服从正态分布，和均未知。现测得16只元件的寿命如下: 159 280 101 212 224 379 179 264 222 362 168 250 149 260 485 170 问是否有理由认为元件的平均寿命大于225小时? 解:未知，按题意作如下假设 H0:0225， H1:0225 取0.05。在Matlab实现 >〉 X=[159 280 101 212 224 379 179 264 222 362 168 250 149 260 485 170］; ％注意： 此处数据X只能为向量而非矩阵 〉> ［h，sig］=ttest（X，225，0。05，1） h = 0 sig = 0.2570 结果表明,h = 0,即在显著水平为0。05的情况下,不能拒绝原假设，认为元件的平均寿命不大于225小时. 8。8.2 两个正态总体均值差的检验（t检验) 函数：ttest2 %具有相同方差的两个正态总体样本均值的比较 格式：[h, sig， ci］ = ttest2 （X, Y) [h， sig, ci] = ttest2 （X， Y， alpha） ［h, sig， ci］ = ttest2 （X， Y, alpha， tail） 说明：原假设为XY 备择假设为：当tail = 0时，表示XY(缺省） 当tail = 1时，表示XY 当tail = —1时，表示XY (X，Y分别表示X，Y的期望) h，sig,ci的含义与前面相同。 例8—38 在平炉上进行一项试验以确定改变操作的建议是否会增加钢的得率，试验是在同一只平炉上进行的。每炼一炉钢时除操作方法外，其它条件都尽可能做到相同。先用标准方法炼一炉，然后用建议的新方法炼一炉，以后交替进行,各炼10炉,其得率分别为： （1)标准方法：78。1 72。4 76.2 74。3 77.4 78.4 76。0 75.5 76.7 77。3 （2）新方法： 79。1 81.0 77.3 79.1 80。0 79。1 79.1 77.3 80。2 82。1 设这两种方法相互独立，且分别来自正态总体N （1,2)和N (2,2)，1、2、2均未知。问建议的新操作方法能否提高得率?（取0.05） 解:（1)建立假设 H0：1=2 H1:1〈2 (2）Matlab实现 〉〉 X=［78。1 72.4 76.2 74.3 77。4 78。4 76。0 75。5 76.7 77。3]； >> Y=[79.1 81.0 77。3 79。1 80.0 79.1 79。1 77。3 80。2 82.1]； 〉> ［h,sig，ci]=ttest2（X，Y，0.05，-1) h = 1 sig = 2。1759e—004 28 ci = —Inf -1.9083 结果h = 1，表明在0.05的显著水平下，可以拒绝原假设，即认为建议的新操作方法较原方法优。 8.8。3 一般性的连续总体的检验 8。8.3.1 两个总体一致性的检验——秩和检验 函数：ranksum 格式：P = ranksum （X, Y， alpha) ［p， h］ = ranksum （X， Y, alpha) 说明：P——两个总体样本X和Y为一致的显著性概率，若P接近于0,则不一致较明 显,可对原假设提出质疑。 X、Y—-两个总体的样本，可以不等长 alpha——给定的显著性水平 h-—检验结果： h = 0--表示X与Y的总体差别不显著 h = 1——表示X与Y的总体差别显著 例8-39 某商店为了确定向公司A或公司B购买某种商品,将A和B公司以往的各次进货的次品率进行比较，数据如下所示.设两样本独立,问两公司的商品的质量有无显著差异。设两公司的商品次品的密度最多只差一个平移，取0.05。 A：7.0 3。5 9。6 8.1 6.2 5.1 10。4 4.0 2。0 10.5 B:5。7 3。2 4。1 11.0 9。7 6.9 3。6 4.8 5.6 8.4 10。1 5.5 12.3 解：分别以A、B记公司A、B的商品次品率总体的均值。该问题是在显著水平0.05下检验假设： H0:A＝B H1：A≠B Matlab实现 >> A=［7.0 3.5 9。6 8.1 6。2 5。1 10。4 4.0 2。0 10.5］； 〉〉 B=[5。7 3。2 4。2 11。0 9。7 6.9 3.6 4.8 5。6 8.4 10。1 5.5 12.3］; 〉〉 ［p,h]=ranksum(A，B,0.05) p = 0.8041 h = 0 结果说明，两样本总体均值相等的概率为0.8041，并不接近于0,且h = 0说明可以接受原假设，即认为两个公司的商品质量无明显差异. 8.8.3。2 中值检验 在假设检验中还有一种检验方法为中值检验，在一般的教学中不一定介绍，但它在实际中也是被广泛应用到的。在Matlab中提供了这种检验的函数，函数的使用方法简单，本文只给出函数介绍。 1. signrank函数 函数：signrank %两个总体中位数相等的假设检验——符号秩检验 格式：P = signrank (X， Y, alpha） [p, h］ = signrank （X, Y, alpha） 说明：P——两个样本X和Y的中位数相等的概率。如果P接近于0，可对原假设质疑。 X、Y-—两个总体的样本，长度必须相同. h-—检验结果: h = 0——表示X与Y的中位数相差不显著； h = 1——表示X与Y的中位数相差显著。 alpha——给定的显著水平。 2. signtest函数 函数:signtest ％两个总体中位数相等的假设检验——符号检验 29 格式:P = signtest (X， Y， alpha) [p， h］ = signtest (X， Y， alpha） 说明：与signrank函数相同。 8.9 方差分析与回归分析 8.9。1 方差分析 在科学试验和生产实践中,影响一事物的因素往往是很多的。例如在化工生产中，有原料成分、原料剂量、催化剂、反应温度、压力、溶液浓度、反应时间、机器设备及操作人员的水平等因素.每一因素的改变都有可能影响产品的数量和质量。有些因素影响较大，有些较小.为了使生产过程得以稳定、保证优质、高产地进行，就有必要找出对产品质量有显著影响的那些因素，为此需进行试验。方差分析就是根据试验的结果进行分析，鉴别各个有关因素对试验结果影响的有效方法。 在试验中,将要考察的指标称为试验指标，影响试验指标的条件称为因素。因素可分为两类：一类是人们可以控制的（可控因素)；一类是人们不能控制的。例如反应温度、原料剂量、溶液浓度等是可以控制的，而测量误差、气象条件等一般是难以控制的。以下所说的因素都是指可控因素。因素所处的状态,称为该因素的水平。 8。8。1.1 单因素方差分析 一项试验如果只有一个因素在改变称为单因素试验,如果多于一个因素在改变称为多因素试验。单因素方差分析的基本问题，是比较和估计多个等方差正态总体的均值。 基本模型为： Xij0jij 其中：Xij——观测值矩阵 0j-—其各列为数据矩阵Xij各列均值的矩阵（0j表示中第j列的所有行） ij-—随机扰动矩阵 在Matlab中，单因素分差分析由函数anova1实现。 anova1是一个4×6表： Source——方差来源；SS-—平方和；df—-自由度；MS——均方值；F-—统计量； Prob>F——P的值;Columns—-因素；Error——误差;Total——总和。 函数:anova1 ％单因素试验的方差分析（anova），比较两组或多组数据的均值，返 回均值相等的概率值。 格式：P = anova1 (X） P = anova1 (X， group) 说明：anova1 (X）—-对样本X中的两列或多列数据进行均衡(即各组数据个数相等）的单因素方差分析，比较各列的均值.其返回值P表示X中各列的均值相等的概率值,如果该值接近于0（即P<），则各列均值不相等显著。anova1 （X) 除了给出P外，还输出一个方差表和一个箱形图，箱形图反映了各组数据的特征。 X——数据矩阵，其各列为各样本值。 anova1 (X, group)——处理非均衡数据（即各组数据个数不相等)， X为数组（向量）,从第1组到第r组数据依次排列； group为与X同长度的数组（向量)，标志X中数据的组别（在与X第i组数据相应的位置处输入整数i(i = 1, 2， …， r））。 例8—40 设有3台机器，用来生产规格相同的铝合金薄板。取样，测量薄板的厚度，精确至‰厘米.得结果如下: 机器1：0.236 0.238 0。248 0。245 0。243 机器2：0。257 0。253 0.255 0。254 0.261 机器3：0。258 0。264 0。259 0。267 0.262 检验各台机器所生产的薄板的厚度有无显著的差异? 30 解：此问题为单因素试验检验问题。 在Matlab中实现为： 〉> X=［0.236 0.238 0。248 0。245 0.243； 0。257 0。253 0。255 0.254 0。261; 0.258 0.264 0。259 0。267 0.262］; >〉 P=anova1（X') P = 1。3431e-005 anova1函数的数据结构结果图形 anova1函数的箱形图 由于P的值接近于0，故拒绝原假设，认为各台机器生产的薄板厚度有显著的差异。 例8—41 用4种工艺生产灯泡，从各种工艺制成的灯泡中抽取了若干个测量其寿命（小时），结果如下表，试推断这几种工艺制成的灯泡寿命是否有显著差异. 工艺 A1 A2 A3 A4 序号 1 1620 1580 1460 1500 2 1670 1600 1540 1550 3 1700 1640 1620 1610 4 1750 1720 1680 5 1800 解：因素A的不同水平下试验数据个数（样本容量）不同，故用函数anova1 (X， group）。 为方便计，将表中数据除以10，再减去150后输入。Matlab实现： >> X=［12 17 20 25 30 8 10 14 22 —4 4 12 0 5 11 18]; >> g=[1 1 1 1 1 2 2 2 2 3 3 3 4 4 4 4]; 〉> anova1(X,g） ans = 31 0.0331 因为0。01〈p = 0.0331<0.05，所以几种工艺制成的灯泡寿命有显著差异。 8.9.1.2 双因素方差分析 基本数学模型 Xijk0ji0ijijk 其中：Xijk——观测值矩阵 ——样本总体均值（常数均值） 0j——列元素为组均值的矩阵(各行总和为0) i0--行元素为组均值的矩阵（各列总和为0） ij——交互作用项（矩阵）(各行、列的总和为0） ijk--随机干扰矩阵 在Matlab中,双因素方差分析用函数anova2实现。 函数:anova2 格式：p = anova2 (X, reps） 说明：anova2 （X， reps)进行平衡的双因素试验的方差分析，用来比较样本X中多行或多列数据的均值，不同列中的数据表示单因素的变化情况,不同行中的数据表示另一因素的变化情况,如果每个行列对（“单元\")有多于一个的观测值，则变量reps指出每一个“单元”观测值的数目，每一单元包含reps行。 例8-42 一火箭使用了4种燃料，3种推进器作射程试验，每种燃料与每种推进器的组合各发射火箭两次,测得结果如下： 推进器B B1 B2 B3 燃料A 56。2000 58。2000 65。3000 A1 52.6000 41。2000 60.8000 51。6000 49.1000 54。1000 A2 42.8000 50.5000 48。4000 60。1000 39.2000 70.9000 A3 73.2000 40。7000 58。3000 A4 75.8000 58.2000 48.7000 32 71.5000 51。0000 考察推进器和燃料这两个因素对射程是否有显著的影响？ 解：在Matlab中应用anova2函数来解决此问题。 >> X=［58.2000 56。2000 65.3000 52。6000 41.2000 60.8000 49。1000 54。1000 51.6000 42。8000 50。5000 48.4000 60。1000 70。9000 39。2000 58。3000 73。2000 40。7000 75。8000 58.2000 48。7000 71。5000 51.0000 41.4000］； 〉〉 p=anova2(X，2） p = 0.0035 0。0260 0.0001 41。4000 anova2函数分析图 由结果可知，各试验均值相等的概率都是小概率，故可拒绝概率相等假设，即认为不同燃料或不同推进器下的射程有显著差异。也就是说，燃料和推进器对射程的影响都是显著的。 8.9.2 回归分析 在客观世界中普遍存在着变量之间的关系。变量之间的关系一般来说可分为确定性的与非确定的两种。确定性关系是指变量之间的关系可以用函数关系来表达的。另一种非确定性的关系即所谓相关关系。例如人的身高与体重之间存在着关系，一般来说,人高一些,体重要重一些，但同样身高的人的体重往往不相同。人的血压与年龄之间也存在着关系，但同年龄的人的血压往往也不相同。气象中的温度与湿度之间的关系也是这样。这是因为涉及的变量（如体重、血压、湿度）是随机变量。上面所说的变量关系是非确定性的。 回归分析是研究相关关系的一种数学方法，是用统计数据寻求变量间关系的近似表达式—-经验公式，它能帮助从一个变量取得的值去估计另一变量所取的值。 在Matlab中,统计回归问题主要由函数polyfit实现。 函数：polyfit 格式：[p, s］ = polyfit （X, Y, n） 说明：p——拟合n次多项式系数向量。 s—-为拟合多项式系数向量的结构信息，返回用函数ployval （ )获得的误差分 析报告。 例8-43 为了研究某一化学反应过程中，温度x对产品得率y的影响，测得数据如下: 温度x 100 110 120 130 140 150 160 170 180 190 得率y 45 51 54 61 66 70 74 78 85 89 试作y = a + bx型的回归. 解：在Matlab中实现 >> x=［100 110 120 130 140 150 160 170 180 190］; 〉〉 y=[45 51 54 61 66 70 74 78 85 89]; >> ［p,s]=polyfit(x,y,1） p = 0。4830 —2。7394 s = R： ［2x2 double］ 33 df: 8 normr： 2。6878 结果说明回归直线为y = -2。7394 +0.4830x 拟合曲线为 >> xi=100:10:190; 〉〉 pi=polyval（p,xi); 〉> plot（x，y，'＊’，xi，pi） 拟合曲线图 34 因篇幅问题不能全部显示，请点此查看更多更全内容