MATLAB在概率统计中的应用
目录
一概率部分
1. 随机变量概率分布的概率计算以及数字特征……………………………………………21.1 随机变量概率分布的概率计算……………………………………………………………21.2 随机变量概率分布的数字特征………………………………………………………… 7
二统计部分
2．数理统计的基础概念…………………………………………………………………………103. 参数估计…………………………………………………………………………………… 114. 假设检验……………………………………………………………………………………165.一元线性回归分析…………………………………………………………………………… 271 随机变量概率分布的概率计算以及数字特征
1.1 随机变量概率分布的概率计算
在MATLAB中列举了多种常见的概率分布，给出了这些概率分布的分布密度函数、分布函数、逆分布函数、随机数发生函数等等，在这一节中，主要研究的是常见概率分布的数字特征（数学期望，方差，协方差以及相关系数）和一些概率的计算
MATLAB中列举的离散型随机变量包括：离散均匀分布、二项分布、泊松分布、几何分布、 超几何分布、负二项分布（Pascal分布）：连续型随机变量包括：连续均匀分布、指数分布、正态分布、对数正态分布、 2分布、非中心 2分布、t分布、非中心t分布、F
χ χ
分布、非中心F分布、β分布、γ分布、Rayleigh分布、Weibull分布。
下表是对这20种分布中的常见分布在Matlab中的应用的总结
表一 常见分布的密度函数在x处的值
分布类型 函数 函数调用格式 备注
名 称 名称
p=normpdf（X，MU， 2
正态分布 normpdf 计算正态分布N(μ,σ )的密度函数在x处
SIGMA） 的值，其中参数SIGMA是μ，MU是σ二项分布 binopdf p=binopdf(x,n,p)
均匀分布 unifpdf p=unifpdf（x，a，b） 计算均匀分布U[a,b]的密度函数在x处的值
几何分布 geopdf p=geopdf(a,p)
超几何分 hygepdf p=hygepdf(x,m,k,n)
布
指数分布 exppdf p=exppdf（x，λ） 计算指数分布的密度函数在x处的值泊松分布 poisspdf p=poisspdf(x,λ)
t分布 tpdf p=tpdf（x，n） 计算t分布的密度函数在x处的
2分布 chi2pdf p=chi2pdf（x，n） 计算 2分布的密度函数在 处的值χ χ x
F分布 fpdf p=fpdf（x，n1，n2） 计算F分布的密度函数在x处的值 表二 运用matlab计算常见分布的分布函数
分布类型 函数 函数调用格式以及意义 备注
名 称 名称
若 2 ,计算 可用
p=normcdf（x， X~N(μ,σ ) P{X >x}
mu，sigma） p1=normcdf(x,mu,sigma)
p=1-p1
若 2 ,计算 可用
X~N(μ,σ ) P{x1
p1=normcdf(x1,mu,sigma)
机变量落在(?∞,x]的概
p2=normcdf(x,mu,sigma)[1]
率，其中mu是参数μ， 2
p=p2-p1[1]
sigma是参数σ
或者p=normspec([x1x2],mu,sigma)
p=binocdf(x,n,p) 若求P(ξ
机变量落在(?∞,x]的概
率 若求P(ξ≥x)则p=1-binocdf(x-1,n,p)
若X ~U[a,b],计算P{X >x}可用
Y=unifcdf（x，a，b p1=unifcdf(x,a,b)
p=1-p1
均匀分布 unifcdf 计算服从均匀分布的随 若X ~U[a,b],计算P{x
1 2
机变量落在(?∞,x]的概 p1=normcdf(x,a,b)
1
率 p2=normcdf(x2,a,b)
p=p2-p1
p=geocdf(a,p) 若求P(ξ
机变量落在(?∞,x]的概
率 若求P(ξ≥x)则p=1-geocdf(k-1,p) p=hygecdf(x,m,k,n) 若求P(ξ
随机变量落在(?∞,x]的
概率 若求P(ξ≥x)则p=1-hygecdf(x,m,k,n)
若X服从参数为λ的指数分布,计算P{X>x}
可用
Y=expcdf（x，λ） p1=unifcdf(x,λ)
计算服从指数分布的随 p=1-p1
指数分布 expcdf 机变量落在 的概 若X服从参数为λ的指数分布,计算
(?∞,x] 可用
P{x1
率 p1=normcdf(x,λ)
1
p2=normcdf(x2,λ)
p=p2-p1
p=poisscdf(x,λ) 若求Pξ(
机变量落在(?∞,x]的概
率 若求P(ξ≥x)则p=1-poisscdf(x,m,k,n)
若X服从自由度为n的t分布,计算P{X>x}
可用
Y=tcdf（x，n） p1=chi2cdf(x,n)
p=1-p1
t分布 tcdf 若X服从自由度为n的t分布,计算
计算服从t分布的随机 可用
P{x1
变量落在(?∞,x]的概率 p1=chi2cdf(x,n)
1
p2=normcdf(x2,n)
p=p2-p1
若X服从自由度为n的 2分布,计算
χ P{X>x}
Y=chi2cdf（x，n） 可用
p1=chi2cdf(x,n)
p=1-p1
χ2分布 chi2cdf 计算服从χ2分布的随
若X服从自由度为n的χ2分布,计算
机变量落在(?∞,x]的概
P{x
率 1 2
p1=chi2cdf(x1,n)
p2=normcdf(x2,n)
p=p2-p1
若X服从自由度为（
n1，n2）的F分布,计算
P{X>x}可用
Y=fcdf（x，n1，n2）
p1=fcdf(x,n)，p=1-p1
若X服从自由度为（n1，n2）的F分布,计算F分布 fcdf 计算服从F 分布的随
P{x1
机变量落在(?∞,x]的概
率 p1=fcdf(x1,n)
p2=normcdf(x2,n)
p=p2-p1
例1：设X~N(3,1.52)
（1） 求X的密度函数在x=2是的值
（2） 求P{x<1},P{1
解：(1) >>p=normpdf(2,3,1.5)
p=
0.2130
所以X的密度函数在x=2是的值是0.2130
（2）令p1=P{x<1}
>>p1=normcdf(1,3,1.5) 结果：p1= 0.0912
令p2=P{1
方法一：
>>p=normcdf(1,3,1.5);
q=normcdf(3,3,1.5);
p2=q-p
结果：p2=0.4088
方法二：>>p2=normspec([1,3],3,1.5)
结果：p2=0.40879
CriticalValue即临界值
Density 即密度
图中蓝色部分表示随机变量X~N(3,1.52)，变量X在[1.5，3]的概率为0.48079由蓝色曲线与横轴围成的部分的概率为1
令p3=P{x?4>2}，因此p3=1?P{x?4≤2}=1?P{2≤x≤6}
>>p3=1-normcdf(6,3,1.5)+normcdf(2,3,1.5)
结果：p3= 0.2752
或者
>>p3=1-normspec([2,6],3,1.5)
结果：p3 =0.2752
例2：生产某种产品的废品率为0.1，抽取20件产品，初步检查已发现有2件废品，问这20件中，废品不少于3件的概率。[4]
解：设抽取20件产品中废品的个数为ξ，则ξ~B(20,0.1)，由于初步检查已发现有2件废品，说明已知20件产品中废品数ξ≥2，因此是求在给定ξ≥2的条件概率
于是P{ξ≥3|ξ≥2}=P{ξ≥3,ξ≥2}÷P{ξ≥2}=P{ξ≥3}÷P{ξ≥2}
令P=P{ξ≥3|ξ≥2}
>>p=(1-binocdf(3,20,0.1))/(1-binocdf(2,20,0.1)) 结果：p=0.4115
例3：某人进行射击试验，假定在300米处向目标射击的命中率为0.02，现独立射击500次，问至少命中3发的概率是多少？[4]
解：将每次射击视为一次试验E，500次射击相当于作500重Bernouli试验E500.用ξ 表示E500 击中目标的次数，依题意， ξ 服从参数为 n=500，p=0.02的二项分布b(x,500 ,0.02)，于是，所求概率为
P(ξ ≥3) =1?P(ξ ≤2)=1?P(ξ =0)?P(ξ =1)?P(ξ =2)
由于n足够大，p足够小，所以可以用泊松分布近似，
p = 0.02,n = 500 ,λ = np =10
在matlab中的实现程序为：
>>p=1-poisspdf(0,10)-poisspdf(1,10)-poisspdf(2,10) 结果：p =0.9972
所以P(ξ≥3)≈0.9972
例4：修理某机器所需时间(单位：小时)服从参数λ=0.5指数分布，试求
（1） 修理时间超过2小时的概率是多少？
（2） 若已经持续修理了9小时，问还需要至少一小时才能修好的概率是多少？[4]解：（1）表示修理时间ξ , ξ 服从参数λ =0.5指数分布,实际是求
p{ξ >2}=1?P{ξ ≤2}
程序如下：
>>p=1-expcdf(2,0.5) 结果： p=0.0183
（3） 由指数分布的无记忆性可知
P(ξ >1+9 |ξ > 9)= P(ξ >1)=1? P(ξ ≤1)
程序如
下：
>>p=1-expcdf(1,0.5) 结果：p =0.1353
例5：设随机变量服从U[0,5]上均匀分布，问方程4x2+4ξx+ξ+2=0,有两个不同的实数根的概率是多少？
解：4x2+4ξx+ξ+2=0,有两个不同的实数根，则
?=16ξ2?16(ξ+2)>0,则ξ 2
P(有两个不同的实数根)=p(ξ2)=1+p(ξ
程序如下：>>p=1+unifcdf(-1,0,5)-unifcdf(2,0,5) 结果：p=0.6000
所以有两个不同实数根的概率是0.6
1.2 随机变量概率分布的数字特征
随机变量的分布函数完整地描述了随机变量的统计特征，但在实际问题中，常常不容易求出分布函数，而许多时候只是需要知道它的某些特征就足够了，例如期望，方差，协方差以及分位数等，在MATLAB6.5的工具箱中提供了求20种分布的数字特征的函数，以下是对本科中常见的10中分布进行归纳总结以及对MATLAB程序的纠正，对于求一般的概率分布期望与方差则用其定义进行求解
表一：常见10种分布的数学期望与方差
分布类型 函数 函数调用格式 备注
名 称 名称
正态分布 normstat [M,V]=normstat(MU,SI 2
（1）M是期望，V是方差σ ，不是标准
GMA)
二项分布 binostat [M,V]=binostat(n,p) 方差σ，指数分布的期望为λ^-1,
方差为λ^-2,在MATLAB中的expstate.m的均匀分布 unifstat [M,V]=unifstat(a,b) 程序中有错误：应将m=mu,v=mu.^2改为
v=mu.^-2;这里的m 表示期望，v表示方几何分布 geostat [M,V]=geostat(p) 差，mu表示参数λ
超几何分 hygestat(m，k，n) （2）如果变量ξ ，i=1,2,3,.......,n相互独布 hygestat i
n n n n
指数分布 expstat [M,V]=expstat(λ) 立则E ξ = Eξ，D ξ = Dξ
∑i ∑ i ∑ i ∑ i
i=1 i=1 i=1 i=1
泊松分布 poisstat p=poisstat(λ) n n
E∏ξi=∏Eξi
t分布 tstat [M,V]=tstat(n) i=1 i=1
（3）如果只是单独求期望或者方差，下面χ2分布 chi2stat [MV]=chi2stat(5)
举例说明：
F分布 fstat p=fstat（n1,n2）
求正态分布的均值：
M=normstat(MU,SIGMA)
求正态分布的方差：
V=normstat(MU,SIGMA)
对于其他分布调用格式类似
例3：计算服从均匀分布U[2,4]的方差与期望
在MATLAB中实现为：
>>[M,V]=unifstat(2,4) 结果：M=3，V=0.3333
例4：计算服从参数λ=0.5的指数分布的方差与期望
>>[M,V]=expstat(0.5) 结果：M=2，V=4
例5：设随机变量(ξ,η)~N(1,1;4,9;r)，求E(ξ?2η+1)
解：由边缘分布可知ξ~N(1,4)，η~N(1,9)
E(ξ?2η+1)=Eξ?2Eη+1
>> M1=normstat(1,2);M2=normstat(1,3);E=M1-2*M2+1 结果E=0
例6：设随机变量 服从区间 1 1 上均匀分布，求 的数学期望
ξ [?2,2]
η=πξ
解：因为ξ服从区间 1 1 上均匀分布，所以πξ服从区间 π π 上均匀分布，，所以
[?2,2] [?2,2]
>>M=unifstat(-pi/2,pi/2) 结果M=0，所以η=πξ的数学期望是0
表二：常见10种分布的p分位数
分布类型 函数 函数调用格式
名 称 名称
正态分布 norminv xp=normminv（p，mu，sigma）
二项分布 binoinv xp=binoinv(p,n,p1)
均匀分布 unifinv xp=normminv（a,b）
几何分布 geoinv xp=geoniv(k,p)k=0,1,2……
超几何分 hygeinv x=hygeinv(p,m,k,n)
布 p
指数分布 expinv xp=expinv（λ）
泊松分布 poissinv xp=poissinv（λ）
t分布 tinv xp=tinv(n)
χ2分布 chi2inv x=chi2inv（n）
p
F分布 fstat xp=finv（n1,n2）
例7：某工厂生产的产品中废品率为0.005，任意取出1000件，计算
（1） 其中至少2件废品的概率：
（2） 其中不超过5件废品的概率
（3） 能以0.9以上的概率保证废品件数不超过多少？
解：用x表示取出的1000件中的废品数，则X近似服从二项分布B(1000,0.005)(1) 所求概率为P{x≥2}=1?P{x<2}=1?P{x≤1}
>>p1=1-binocdf(1,1000,0.005) 结果p1=0.9599
所以其中至少2件废品的概率是0.9559
(2) 所求概率为P{x≤5}
>>p2=binocdf(5,1000,0.005) 结果p2=0.6160
(3)由题意，求n，使P{X<=n}>=0.9,我们只求n满足P{X<=n}=0.9即可
>>n=binoinv(0.9,1000,0.005)结果n=8
所以能以0.9以上的概率保证废品数不超过8件
2 数理统计的基础概念
2.1 总体与样本
在统计学中把研究的对象全体称为全体，确切地说，应该是这些对象的某个或某些指标值全体称为总体，为了了解总体的分布，从总体中抽取n个个体，记其指标值为x1,x2,......, xn，称x1,x2,......, xn为总体的一个样本，n为样本的大小或容量
利用MATLAB数理统计工具箱中的函数可以直接计算样本的均值，方差以及其他数字特征
表一：MATLAB中计算样本的数字特征的函数
数字特征 调用函数 备注
若X为向量，返回X中各元素的平均值
样本平均值 M=mean(X) 若X为矩阵，返回X中各列元素的平均值构成的向
量
若X为向量，返回X中各元素的平均值
中位数 N=median(X) 若X为矩阵，返回X中各列元素的平均值构成的向
量
1 n
var(X)= (x ?X)2
∑ i
n?1i=1
样本方差 D=var(X) 若X为向量，则返回向量的样本方差。
若X为矩阵，则D为X的列向量的样本方差构成的
行向量
返回向量（矩阵）X的样本标准差（置前因子为 1 ）样本标准方差 D=std(X) n?1 1 n
即：std= x ?X
n?1∑ i
i=1
协方差 C=cov(X) 返回向
量X的协方差
相关系数 返回列向量X,Y的相关系数，等同于corrcoef([X Y])
P=corrcoef(X,Y)
hist(x,k) 众数通常是从频数直方图上近似看出来，因此求众数
hist(x,nb) 调用的是画直方图的函数，其中x是由数据组成的向
众数 [n,x]=hist(x,k) 量，参数k表示将区间k等分，默认值为10，参数
[n,x]=hist(x,nb) nb为事先给出的小区间的区间图，命令1和2能画直
方图，命令3和4不画出直方图，其输出参数n为落
入每个小区间的频数，x为区间的中点
例子8：统计某班30人某门课程某次考试成绩的分布情况，考试成绩为
9375839391858482777695948991888683968179
9778756769688483817566857094838280787473
767086769089716686738084797877635355
求平均成绩与成绩的中位数以及方差和标准差和众数[5]
解：>>x=[9375839391858482777695948991888683968179...
9778756769688483817566857094838280787473...
767086769089716686738084797877635355];
M=mean(x)
N=median(x)
D1=var(x)
D2=std(x)
hist(x,4)
结果：
M=79.9138
N=80.5000
D1=
94.0100
D2=9.6959
所以平均成绩是79.9138，成绩中位数是80.5000，方差是94.0100，标准差是9.6959从直方图中可以近似看出众数为80
3 参数估计
3.1 矩估计
设总体X具有分布函数F(x;θ)，θ为p维待估参数，X的k阶中心矩为
k k
，样本 的k阶中心矩为 1 n ，样本的k阶中心矩m =E(X?E(X)) X,X ,...,X
k 1 2 n A = (X?X)
k n∑ i
i=1
的计算可用函数moment，其调用格式是B=moment(X,k)
当p=1时，由E(X)=X求出θ的矩估计，其中样本均值X可用函数mean(X)计算当p≠1时，令
? E(X)=X
? k=2,3,....p
?mk(θ)=Ak,
解此方程组，即可得到θ的矩估计
?1
, 0≤x≤θ
?θ
例9：设变量X服从区间[0,θ]上的均匀分布，即分布函数为p(x;θ)=?
?
?0, 其他
?
现得样本值为1.3，0.6，1.7，2.2，0.3，1.1，试用矩估计法求参数θ的估计
θ 1
解：X=E(X)=∫x dx
0 θ
>>Symsx;
Symssita;
θ 1
y=int(x./sita,0,sita) （求 x dx的值，其中记θ为sita）
∫ θ
0
θ 1 1
结果：y =1/2*sita,所以∫x dx= θ
0 θ 2
>>X=[1.30.61.72.20.31.1];
mean(X)
结果：ans=1.2000
1 ?
所以X=1.2,所以 2θ=1.2，所以θ=2.4
例10：已知某种灯泡的寿命(单位：h)服从正态分布，在某周所生产的该种灯泡中随机抽取10只，测得其寿命为1067，919，1196，785，1126，936，918，1156，920，948.设总体参数为都未知，试用矩估计求出参数μ,σ的估计值
解：? E(X)=X
? 2
E(X?E(X)) =DX, k=2,3,....p
?
? μ=X
2 k=2,3,....p
?σ =DX,
>>X=[1067 919 1196 78511269369181156 920948];
mean(X) %求样本均值
var(X) %求样本方差
结果：
ans =997.1000
ans =1.7305e+004
>>sqrt(1.7305e+004) %求样本均方差
结果
ans =131.5485
所以μ?=997.1000，σ=131.5485?
3.2 极大似然估计与区间估计
在MATLAB中，这两种估计可同时给出,对于二项分布，泊松分布，正态分布，β分布，γ分布，均匀分布，指数分布，威布尔分布在MATLAB中可将这些分布的分布命令字符和函数命令字符fit连接起来(除均匀分布连接后为unifit)，得到这些分布中参数的极大似然估计值以及置信水平为1?α的置信区间.
表一：MATLAB中求极大似然估计与区间估计的函数调用格式
分布 函数功能以及参数的解释名称 函数及其调用格式
参数p的MLE以及置信区间
二项分布 [phat,pci]=binofit(x,n,alpha) 输出参数phat和pci分别表示p的MLE，
以及置信区间
参数λ的MLE以及置信区间，输出参数泊松分布 [lambdahat,lambdaci]=poissfit(x,alpha) lambdahat和lambdaci分别表示
lambda（即λ）的MLE，以及置信区间
（1）参数μ,σ的MLE以及置信区间，
输出参数muhat和muci分别表示μ的
α MLE以及置信区间，sigmahat和sigmaci
[muhat,sigmahat,muci,sigmaci]
未知 分别表示σMLE以及置信区间
=normfit(x,alpha)
正态分布 （2）在使用中发现MATLAB6.5的
normfit.m程序中求σ的极大似然估计
公式用错了，它的原句是
singmahat=std(x)
[muci]= 更正方法为在此M文件的最后一定要加
α normfit1(x,sigma,alpha)
已知 上语句：
singmahat=sqrt(moment(x,2))[2]
参数a和b的MLE以及置信区间，输均匀分布 出参数ahat和aci分别表示a的MLE
[ahat,bhat,aci,bci]=unifit(x,alpha) 以及置信区间，bhat和bci分别表示b
的MLE以及置信区间
参数θ的MLE以及置信区间指数分布 [muhat,muci]=expfit(x,alpha) 输出参数muhat和muci分别表示mu
即μ的MLE以及置信区间
γ分布 [phat,pci]=gamfit(x,alpha,options) 参数a和b的MLE以及置信区间，输 出参数phat和pci分别表示p的MLE，
以及置信区间
威布尔分 参数a和b的MLE以及置信区间，输
布 [phat,pci]=weibfit(data,alpha) 出参数phat和pci分别表示p的MLE，
以及置信区间
注：normfit1的程序编写
function[muci]=normfit1(x,sigma,alpha)
ifnargin<3
alpha=0.05;
end
[m,n]=size(x);
ifmin(m,n)==1
x=x(:);
m=max(m,n);
n=1;
end
muhat=mean(x);
muci=zeros(2,n);
normcrit=norminv(alpha/2);
muci=[(muhat-normcrit*sigma/sqrt(m))(muhat+normcrit*sigma/sqrt(m))];
例11：通常，引力常数的测定
值服从均值为μ、标准差为σ的正态分布，某人在实验中使用金球测定引力常数，6次测定观察值：6.683，6.681，6.676，6.678，6.679，6.672试用极大似然估计法对未知参数μ和σ作出估计并求其置信水平为0.95的置信区间解：用binofit函数进行计算
>>x=[6.683 6.6816.676 6.6786.679 6.672];
[muhat,sigmahat,muci,sigmaci]= normfit(x)
结果：
muhat =6.6782，sigmahat=0.0035，muci= 6.67416.6822，sigmaci =0.0024 0.0095所以μ的极大似然估计值为6.6782，置信区间是[6.6741,6.6822], σ的极大似然估计值为0.0035,置信区间是[0.0024,0.0095]
对常见的20种分布，也可以用函数mle求参数的极大似然估计值以及1-α的置信区间，命令是：
[phat,pci]=mle(‘dist’,x,alpha)
函数中输入’dist’为函数分布的命令字符如norm即表示正态分布，输入参数x的样本数据向量，输入参数alpha的默认值为0.05，注意对于二项分布，命令为
[phat,pci]=mle(‘dist’,x,alpha,n)其中n表示试验次数
对于上述例题的数据，若用mle函数计算
>>x=[6.683 6.6816.676 6.6786.679 6.672]；
[phat,pci]=mle('norm',x,0.05)
结果：phat= 6.6782 0.0035 pci= 6.6751 0.0008
6.6812 0.0063
所以μ的极大似然估计值为6.6782，置信区间是[6.6751,6.6812],σ的极大似然估计值为0.0035,置信区间是[0.0008,0.0063]
例12：随机地从一批钉子中抽取16枚，测得其长度(单位：cm)为2.14，2.10，2.13，2.15，2.13，2.10，2.15，2.12，2.14，2.10，2.13，2.11，2.14，2.11.设钉长分布为正态的，试求总体均值μ的90%置信区间：
（1）若σ未知
（2）若已知σ=0.01cm[4]
解：（1）>>x=[2.142.102.132.152.132.122.132.102.152.122.142.102.132.11
2.142.11];
[muhat, sigmahat,muci,sigmaci]= normfit(x,0.1)
结果：
muhat=2.1250，sigmahat =0.0166，muci=2.1175 2.1325，sigmaci =0.0133 0.0246所以总体均值μ的90%置信区间是[2.1175,2.1325]
(2)解：>>x=[2.142.102.132.152.132.122.132.102.152.122.142.102.132.11
2.142.11];
[muci]= normfit1(x,0.01,0.1)
结果：muci=2.1291 2.1209
所以若已知σ=0.01cm，总体均值μ的90%置信区间是[2.1291,2.1209]
4 假设检验
4.1 正态分布的参数假设检验
正态分布的参数检验包括σ2已知时，单个正态总体的均值μ的假设检验（U检验法），σ2未知，单个正态总体的均值μ的假设检验(t检验法)，σ ,σ 未知，但σ =σ 时，
1 2 1 2两个正态总体均值的估计，μ未知，单个正态总体方差的检验，均值未知时，两个正态总体方差的假设检验，在MTLAB中只提供了前面三种的假设检验，但通过自己编程实现了后面两种情况的检验
表一：单个正态分布的参数假设检验在MTLAB中的实现
检验 提前 函数调用格式 备注
参数 条件
(1)原假设：H0: μ=μ0=m
[h,p,ci,zval]=ztest(x,m,sigma,alpha,tail)
2 p为观察值的概率，当p为小概率时则 (2)若h=0，表示在显著性水平
σ
对原假设提出质疑，ci为真正均值μ的 alpha下，不能拒绝原假设；
已知 1-alpha置信区间，zval为统计量的值，
若h=1，表示在显著性水平alpha
sigma为已知的标准差， 下，可以拒绝原假设。均值
μ [h,p,ci,stats]=ttest(x,m,alpha,tail)
(3)若tail=0，表示备择假设：
p为观察值的概率，当p为小概率时则 H : μ≠μ =m（双边检验）；
2 1 0
σ 对原假设提出质疑，ci为真正均值μ的
1-alpha置信区间，stats为“tstat”的 若tail=1，表示备择假设：
未知 检验统计量T的值及名为“df”的T统
H1: μ>μ0=m（单边检验）；
计量所服从的t分布的自 由度
若tail=-1，表示备择假设：
[h,c,c1]=chitest(x,b,a,tail) H : μ<μ =m（单边检验）。
1 0
方差 μ
σ x为样本数据，a为显著性水平,b为正
未知 态分布的标准方差
注：μ未知，单个正态总体方差的检验[3]
function[h,c,c1]=chitest(X,b,a,tail)
%单个个正态总体方差的检验(卡方检验）
%调用格式：[h,c,c1]=chitest(X,b,a,tail)
%X样本，a为显著性水平,b为正态分布的标准方差
%TAILmustbe''both'',''right'',or''left'',or0,1,or-1.'
%《应用数理统计》（叶慈南编）p119.
ifnargin<3tail=0;end
ifnargin<2a=0.05;end
n=length(X);
c=(n-1)*var(X)/(b^2);
iftail ==0%two-tailedtest
c1=chi2inv(a/2,n-1);c2=chi2inv(1-a/2,n-1);
if(c<=c1)|(c<=c2)
h=1;
else
h=0;
end
elseiftail== 1%rightone-tailedtest
c1=chi2inv(1-a,n-1);
ifc>=c1
h=1;
else
h=0;
end
elseiftail== -1%leftone-tailedtest
c1=chi2inv(a,n-1);
ifc<=c1
h=1;
else
h=0;
end
else
error('stats:ftest:BadTail',...
'TAILmustbe''both'',''right'',or''left'',or0,1,or-1.');
End
例13. 某车间用一台包装机包装葡萄糖，包得的袋装糖重是一个随机变量，它服从正态分布。当机器正常时，其均值为0.5公斤，标准差为0.015。某日开工后检验包装机是否正常，随机地抽取所包装的糖9袋，称得净重为（公斤）
0.497, 0.506, 0.518, 0.524, 0.498, 0.511, 0.52, 0.515, 0.512问机器是否正常？
解：总体μ和σ 已知，该问题是当σ2为已知时，在水平α=0.05下，根据样本值判断μ=0.5还是μ≠0.5。为此提出假设：
原假设： H0: μ=μ0=0.5
备择假设：H1: μ≠0.5
>>X=[0.497,0.506,0.518,0.524,0.498,0.511,0.52,0.515,0.512];
>>[h,sig,ci,zval]=ztest(X,0.5,0.015,0.05,0)
结果显示为：h=1，sig=0.0248 ci=0.50140.5210， zval =2.2444
结果表明：h=1，说明在水平α=0.05下，可拒绝原假设，即认为包装机工作不正常。例14. 某种电子元件的寿命X（以小时计）服从正态分布，μ、σ2均未知。现测得16只元件的寿命如下
159 280 101 212 224 379 179 264 222 362 168 250
149 260 485 170
问是否有理由认为元件的平均寿命大于225（小时）？
解：未知 2，在水平 下检验假设： ：μ<μ =225， ：μ>225
σ α=0.05 H0 0 H1
>>X=[159280101212224379179264222362168250149260485170];
>>[h,sig,ci]=ttest(X,225,0.05,1)
结果显示为：h=0，sig=0.2570 ci=198.2321 Inf
结果表明：H=0表示在水平α=0.05下应该接受原假设H0，即认为元件的平均寿命不大于225小时。
例题15：在正常生产情况下，印花棉布布幅的宽度服从标准方差为0.0048的正态分布，某日选取该种棉布5匹，测得布幅宽度为1.32，1.55，1.36，1.40，1.44，问该日印花棉布布幅宽度的标准差是否正常？（取α=0.05）
解：H0：σ=0.0048， H1：σ≠0.0048
x=[1.32 1.551.361.401.44];
[h,c,c1]=chitest(X,0.0048,0.05,0)
结果：h=0，c=6.7597e+009，c1=2.7004
所以，在α=0.05下，接受原假设，该日印花棉布布幅宽度的标准差正常
表二：双总体正态分布的方差检验在MTLAB中的实现
检验 提前
参数 条件 函数调用格式 备注
[h,p,ci,stats]=ttest2(X,Y,alpha,tail) 原假设： ， (
H0: μ1=μ2 μ1
σ1,σ2未知，
当p为小概率时则对原假设提出 为X为期望值，μ 为Y的期望
质疑，ci为真正均值μ的1-alpha 2
均值 但σ1=σ2
μ 置信区间。stats为“tstat”的检验 值)
统计量T的值及名为“df”的T统
计量所服从的t分布的自由度 若h=0，表示在显著性水平
alpha下，不能拒绝原假设；
若h=1，表示在显著性水平
alpha下，可以拒绝原假设。
若tail=0，表示备择假设：
[h,f,f1]=ftest0(X,Y,a,tail) H : μ ≠μ （双边检验）；
μ 1 1 2
方差
未知 若tail=1，表示备择假设：σ X,Y为样本，a为显著性水平
H : μ >μ （单边检验）；
1 1 2
若tail=-1，表示备择假设：
H : μ <μ （单边检验）。
1 1 2
注：函数ftest0的程序编写如下：
function[h,f,f1]=ftest0(X,Y,a,tail)
%两个正态总体方差的检验(F检验）
%调用格式：[h,f,f1]=ftest0(X,Y,a,tail) ：《应用数理统计》（叶慈南编）p122.
%X,Y为样本，a为显著性水平
ifnargin<4tail=0;end
ifnargin<3a=0.05;end
n1=length(X);n2=length(Y);
f=var(X)/var(Y);
iftail==0%two-tailedtest
f1=finv(1-a/2,n1-1,n2-1);f2=finv(a/2,n1
-1,n2-1);
if(f1<=f)|(f<=f2)
h=1;
else
h=0;
end
elseiftail==1%rightone-tailedtest
f1=finv(1-a,n1-1,n2-1);
iff<=f1
h=1;
else
h=0;
end
elseiftail==-1%leftone-tailedtest
f1=finv(a,n1-1,n2-1);
iff<=f1
h=1;
else
h=0;
end
else
error('stats:ftest:BadTail',...
'TAILmustbe''both'',''right'',or''left'',or0,1,or-1.');
end
例16. 在平炉上进行一项试验以确定改变操作方法的建议是否会增加钢的产率，试验是在同一只平炉上进行的。每炼一炉钢时除操作方法外，其他条件都尽可能做到相同。先用标准方法炼一炉，然后用建议的新方法炼一炉，以后交替进行，各炼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(μ,σ)和N(μ2 ,σ)，μ、μ2 、σ2均未知。
1 2 1 2
问建议的新操作方法能否提高产率？（取α=0.05）
解：先做方差相等的检验H0：σ =σ ，H1：σ ≠σ
1 2 1 2
>>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,f,f1]=ftest0(X,Y,0.05,0)
结果：h=0，f=1.4945，f1=4.0260
因为h=0，所以接受原假设σ1=σ2
在水平α=0.05下检验假设：H ：μ =μ ，H：μ <μ
0 1 2 1 1 2
>>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,p,ci,stats]=ttest2(X,Y,0.05,-1)
结果：
h=1，p=2.1759e-004，ci=-Inf -1.9083，stats= tstat:-4.295 df:18
结果表明：H=1表示在水平α=0.05下，应该拒绝原假设，即认为建议的新操作方法提高了产率，因此，比原方法好。
4.3 分布的检验
随机变量的分布完全描述了随机现象的规律，因此，确定随机现象的分布十分重要，对于分布的检验，主要讨论两大类，第一类是理论分布完全已知且只取有限个值的情形，第二类是理论分布含有未知参数且只取有限个值的情形，这里讨论的理论分布是泊松分布以及正态，对于正态分布的检验，可以运用matlab中kstest.m函数，调用格式为：
[H,pValue,KSstatistic,criticalValue]=kstest(x,CDF,alpha,tail)
说明：（1）原假设为具有相同连续分布.测试结果为H,若H=0,表示应接受原假设;若H=1,表示可以拒绝原假设
（2）P为假设成立的概率,KSstatistic为测试统计量的值，tail=1
对于第一类的情形，主要采用χ2分布，利用皮尔逊定理，自己编出检验的程序程序如下：function[s1,s2]=possion(a,b,c,d,n,alpha)
t1=dot(a,b)./n;
r=size(c);
t2=[poisspdf(0,t1) poisspdf(1,t1) poisspdf(2,t1) poisspdf(3,t1) poisspdf(4,t1) poisspdf(5,t1)1-poisscdf(5,t1)];
s1=sum((d-n.*t2).^2./(n.*t2));
s2=chi2inv(1-alpha,r-1-1);
ifs1>s2
disp('该数据源不服从泊松分布。')
else
disp('该数据源服从泊松分布。')
end
说明：（1）程序还不都完善，程序中的t2需要根据不同的数据进行调整
（2）在计算数据是否服从已知分布时可以先用作图函数cdfplot将数据的分布图画出.例17：根据近63年的观测资料，上海每年夏季（5月至9月）发生暴雨的天数如下：
暴雨天数 0 1 2 3 4 5 6 7 8 ≥9
年份数 4 8 14 19 10 4 2 1 1 0试问上海每年夏季发生暴雨的天数是否服从泊松分布(α=0.05)[4]
解：>>A=[00001111111122222222222222333333333333333333…3444444444455556678];
cdfplot(A)
>>a=[0123456789];
b=[4814191042110];
c=[0123456];
d=[4814191044];
n=63;
alpha=0.05;
[s1,s2]=possion(a,b,c,d,n,alpha)
结果：
该数据源服从泊松分布。s1=2.9090 s2=NaN 11.0705
所以上海每年夏季发生暴雨的天数是否服从泊松分布
例18：某厂生产灯泡，检验光通量（单位：流明）是否服从正态分布(α=0.05)
今测得n=120个灯泡，数据如下：
216203197208206209206208202203
206213218207208202194203213211
193213208208204206204206208209
213203206207196201208207213208
210208211211214220211203216224
211209218214219211208221211218
218190219211208199214207207214
206217214201212213211212216206
210216204221208209214214199204
211201216211209208209202211207
202205206216206213206207200198
200202203208216206222213209219
解：>>A=[216203197208206209206208202203…]（把A的数据以向量的形式输入）cdfplot(A)
[h,p,k,c]=kstest(A,[],0.05,1)
结果：
h=0 p=1 k=0 c=0.1103
所以检验光通量（单位：流明）服从正态分布
4.2 独立性检验
在MATLAB中没有实现独立性检验的程序，通过用列联表的方法进行独立性检验，编程如下：function[chi2,chi21,h,ni,nj]=indenpendencetest0(nij,alpha)
%调用格式：[chi2,chi21,h]=indenpendencetest0(nij,alpha)
%Author: Ji Lin,2007
ifnargin<2
alpha=0.05;
end
[r,s]=size(nij);
n=sum(sum(nij));
ni=sum(nij');
nj=sum(nij);
chi2=n*sum(sum(((nij-ni'*nj/n).^2)./(ni'*nj)));
chi21=chi2inv(1-alpha,(r-1)*(s-1));
ifchi2>=chi21
h='X与Y有关';
else
h='X与Y独立';
end
函数调用格式：[chi2,chi21,h]=indenpendencetest0(nij,alpha)
例19：对150人进行生活满意程度与其婚姻状况的调查，结果如下：
婚 已婚 丧偶
姻 状
满意 况
程度
很满意 22 10
较满意 19 22
一般 11 36
不满意 8 22
问二者是否有关系（α=0.05）?[4]
解：H0：生活满意程度与婚姻状况无关 H1：生活满意程度与婚姻状况有关
>>nij=[2210;1922;1136;822];
[chi2,chi21,h,ni,nj]=indenpendencetest0(nij,0.05)
结果：
chi2=19.3237 chi2=7.8147 h=X与Y有关 ni=32 41 47 30
nj=60 90
所以可以知道生活满意程度与婚姻状况有关
例18：3个工厂生产同一只种产品，产品分1，2，3三个等级，从3个厂中分别随机抽取了109件、100件、91件产品，等级状况如下：
工
厂 甲 乙 丙
等
级
1 58 38 32
2 28 44 45
3 23 18 14
问各工厂产品质量是否一致（α=0.05）[4]
解：各工厂产品质量一致，可以认为是假设
H0：工厂与质量等级无关 H1：工厂与质量等级有关
>>nij=[583832;284445;231814];
[chi2,chi21,h,ni,nj]=indenpendencetest0(nij,0.05)
结果：
chi2=13.5862 chi21=9.4877 h=X与Y有关 ni=128 117 55
nj=109 100 91
5 回归分析
5.1 一元线性回归分析
一元线性回归分析的回归模型为y=b +b x+ε，其中b ，b为未知参数，b 称为常
0 1 0 1 0数项， 称为回归系数， 称为随机误差，并假定 2
b1 ε ε ~N(0,σ )
在MATLAB中有几种函数可供选择进行回归分析，总结如下表：
函数名 调用格式 作用 备注
polyfit a=polyfit(x,y,1) ?
可求出回归系数b ，b a(1)= ,a(2)= ?
b b
0 1 1 0polyval y0=polyval(a,x0) 计算给定x0时的预测值y0
Polytool h= 它产生的是交互式画面，给出拟合
polytool(x,y,1,alp 直线 ，并用两条红线描出
y=b0+b1x
ha)
1?α的置信区间
regress x0=ones(n,1); 得到回归分析中用到的所有信息 stats表示检验回归模
[b,bint,r,rint,stats] b表示回归系数的最小二乘估计， 型的统计量，包括3个
= ? ? 数值；第1个是复相关
b(1)= b ,b(2)= b
regress(y＇，[x 1 0 系数，第2个是检验回
0
x＇]，alpha) bint表示回归系数的1?α的置信 归方差显著性的F统计
区间 量的值，第3个是对应
r表示残差向量 于所得F统计量的概率
rint表示残差向量置信区间 P,当p≤α时，线性方
stats表示检验回归模型的统计量 程有显著意义，当
p>α，线性方程无显
著意义
例20：炼铝厂测得吕的硬度x与抗张强度y的数据如下：
X 68 53 70 84 60 72 51 83 70 64Y 288 298 349 343 290 354 283 324 340 286（1） 求y对x的回归方程
（2） 检验回归方程的显著
性（α=0.05）；
（3） 求y在x=65处的预测区间（置信度为0.95）[4]
解：（1）>>x=[68 53 708460 72 51837064];
y=[288298349343290354283324340286];
a=polyfit(x,y,1)
结果：
a=1.8007 193.9508
所以y对x的回归方程为y=193.9508+1.8007x
利用函数regress的程序和结果如下：
>> x0=ones(10,1);
[b,bint,r,rint,stats]=regress(y',[x0,x'],0.05)
结果
b=193.9508 1.8007 bint =86.0395 301.8620
0.2209 3.3805
r=
-28.4004 8.6106
28.9982 -2.2120
-11.9945 30.3967
-2.7880 -19.4113
19.9982 -23.1974
rint =
-75.3185 18.5178
-37.8517 55.0729
-17.4575 75.4539
-47.3176 42.8935
-62.4428 38.4537
-14.9365 75.7299
-47.8693 42.2934
-62.3445 23.5219
-29.9330 69.9293
-71.8906 25.4957
stats= 0.4634 6.9091 0.0302
由于b1的置信区间（0.2209 3.3805）不包含0，所以认为线性回归方程显著（3）>> x=[6853 70 846072 51 837064];
y=[288298349343290354283324340286];
a=polytool(x,y,1,0.05)
结果：
a=
310.0002311.0002 312.0002313.0002
从图中可以得到y在x=65处的预测区间是(310.9982±55.0075)
5.2 多元线性回归分析
很多情况下要用多元线性回归的方法才能更好地描述变量间的关系，只是多元线性回归的方法复杂写，计算量很大，因此利用matlab软件很解决了这一问题，对于多元线性回归主要讨论两种情况，如下表
求回归方程的 回归方程显著 计算在x的 备注
多元线性回归 函数调用格式 性检验的函数 预测值的函
调用格式 数调用格式 [b,bint,r,rint,stats] 为了大概确
=regress(y,x,alpha) 定n的值，可若因变量y与 p个自变量x Y=ployval(a,x
i 其中y为列向量，x为
n×(p+1) ) 以通过作图
之间的关系是 维矩阵其第一列必须全为1，b为 或者 函 数 scatter
回归系数，bint为回归系数的1?α 画出散点图y=β +βx +...+β x +ε ploytool(x,y,n
0 1 1 p p 的置信区间，若置信区间包括零， scatter=(x,y,’r’
认为x对y的影响不显著，不包括 ,alpha),为了 );
零则显著 大概确定n的 gridon若因变量y与 是n次多项式 值，可以通过
x a=plotfit(x,y,n)
关系 参数x，y为要拟 作 图 函 数
n scatter画出散
y=β +β x+...+β x +ε 合的数据，n为拟
0 1 n [b,bint,r,rint,
合次数，输出参数 stats] 点图
令 2 n
x1=x,x2 =x ,...xn =x scatter=(x,y,’r’
a是回归系数 =regress(y,x );
则转化为 ,alpha)
gridon
y=β0+β1x1+...+βpxp+ε
例21：下表列出了某地区家庭人均鸡肉年消费量Y与家庭月均收入X，鸡肉价格P1（元/千克
），猪肉价格P（元/千克）与牛肉价格P（元/千克）的相关数据。
2 3
表1 某地区家庭消费肉类情况表
年份 Y/千克 X/元 P P P
1 2 3
1980 2.78 397 4.22 5.07 7.83
1981 2.99 413 3.81 5.2 7.92
1982 2.98 439 4.03 5.4 7.92
1983 3.08 459 3.95 5.53 7.92
1984 3.12 492 3.73 5.47 7.74
1985 3.33 528 3.81 6.37 8.02
1986 3.56 560 3.93 6.98 8.04
1987 3.64 624 3.78 6.59 8.39
1988 3.67 666 3.84 6.45 8.55
1989 3.84 717 4.01 7 9.37
1990 4.04 768 3.86 7.32 10.61
1991 4.03 843 3.98 6.78 10.48
1992 4.18 911 3.97 7.91 11.4
1993 4.04 931 5.21 9.54 12.41
1994 4.07 1021 4.89 9.42 12.76
1995 4.01 1165 5.83 12.35 14.29
1996 4.27 1349 5.79 12.99 14.36
1997 4.41 1449 5.67 11.76 11.76
1998 4.67 1575 6.37 13.09 13.09
1999 5.06 1759 6.16 12.98 12.98
2000 5.01 1994 5.89 12.8 12.82001 5.17 2258 6.64 14.1 14.12002 5.29 2478 7.04 16.82 16.82（1）求回归方程
（2）检验回归方程的显著性
解：>>y=[2.78
2.99
2.98
……
5.29];
x=[13974.225.077.83
14133.815.27.92
14394.035.47.92
……..
124787.0416.8216.82];
[b,bint,r,rint,stats]=regress(y,x,0.05)
b=
3.6564
0.0011
-0.5169
0.1229
0.0463
bint=
2.7145 4.5983
0.0006 0.0016
-0.8533 -0.1805
-0.0547 0.3006
-0.0608 0.1533
r=
-0.1336
-0.1740
-0.1245
-0.1047
-0.1487
-0.1137
0.0659
0.0270
0.0499
0.1440
0.1116
0.1504
0.0361
0.2672
0.0276
-0.1419
-0.1945
0.0408
0.2938
0.3838
-0.0435
-0.0172
-0.4017
rint=
-0.4701 0.2028
-0.5655 0.2175
-0.5106 0.2615
-0.5021 0.2927
-0.5447 0.2473
-0.4992 0.2718
-0.2928 0.4245
-0.3639 0.4179
-0.3529 0.4528
-0.2595 0.5475
-0.2659 0.4891
-0.2103 0.5110
-0.3246 0.3968
-0.0821 0.6165
-0.3526 0.4078
-0.4754 0.1917
-0.5389 0.1499
-0.3511 0.4327
-0.0497 0.6373
0.0400 0.7276
-0.4079 0.3208
-0.3536 0.3193
-0.6915 -0.1120
stats=
0.9412 72.0222 0.0000
由于β的0.95的置信区间为（0.0006 0.0016）不包含0，故拒绝H01：β=0，
1 1由于β 的0.95的置信区间为（-0.8533 -0.1805）不包含0，故拒绝H02：β =0，
2 2由于β 的0.95的置信区间为（-0.0547 0.3006）包含0，故接受H03
：β =0，
3 3由于β 的0.95的置信区间为（-0.0608 0.1533）包含0，故接受H04：β =0
4 4所以P2，P3对Y的影响不显著，X，P1对Y的影响显著
所以回归方程：Y =3.6564+0.0011X?0.5169P1