Matlab统计工具箱的应用

Matlab统计工具箱的应用

一、统计的基本命令

1、数据的录入、保存和调用

例1上海市区社会商品零售总额和全民所有制职工工资总额的数据如下

方法1（将单个变量存到数据文件.mat中）

（1）年份数据以1为增量，用产生向量的方法输入。

命令格式： x=a:h:b t=78:87

（2）分别以x和y代表变量职工工资总额和商品零售总额。

x=[23.8,27.6,31.6,32.4,33.7,34.9,43.2,52.8,63.8,73.4] y=[41.4,51.8,61.7,67.9,68.7,77.5,95.9,137.4,155.0,175.0]

（3）将变量t、x、y的数据保存在文件data中: save data t x y

（4）进行统计分析时，调用数据文件data中的数据: load data

方法2(将数据矩阵存入数据文件中)

（1）输入矩阵：

data=[78,79,80,81,82,83,84,85,86,87,88;

23.8,27.6,31.6,32.4,33.7,34.9,43.2,52.8,63.8,73.4;

41.4,51.8,61.7,67.9,68.7,77.5,95.9,137.4,155.0,175.0]

（2）将矩阵data的数据保存在文件data1中：save data1 data

（3）进行统计分析时，先用命令： load data1

调用数据文件data1中的数据，再用以下命令分别将矩阵data的第一、二、三行的数据赋给变量t、x、y： t=data(1,:); x=data(2,:); y=data(3,:);

若要调用矩阵data的第j列的数据，可用命令：data(:,j)

Chengxu2；chengxu22；

方法3(直接将数据存入文本文件中，然后调用)

（1）建立data.txt文件报保存；

（2）load data.txt

2、基本统计量

对随机变量x ，计算其基本统计量的命令如下：

均值：mean(x) 中位数：median(x) 标准差：std(x) 方差：var(x) 偏度：skewness(x) 峰度：kurtosis(x) 协方差：cov （x,y ） 相关系数：corrcoef （x,y ） 3、常见概率分布的函数 常见的几种分布的命令字符为：

正态分布：norm 指数分布：exp 帕松分布：poiss β分布：beta 威布尔分布：weib 2χ分布：chi2 t 分布：t F 分布：F Matlab 工具箱对每一种分布都提供五类函数，其命令字符为：

概率密度：pdf 概率分布：cdf 逆概率分布：inv 均值与方差：stat 随机数生成：rnd

（当需要一种分布的某一类函数时，将以上所列的分布命令字符与函数命令字符接起来，并输入自变量（可以是标量、数组或矩阵）和参数即可. 如对均值为mu 、标准差为sigma 的正态分布，举例如下： 例 画出正态分布)1,0(N 和)2,0(2N 的概率密度函数图形.

在Matlab 中输入以下命令：

1)密度函数：p=normpdf (x,mu,sigma ) (当mu=0,sigma=1时可缺省) x=-6:0.01:6; y=normpdf(x); z=normpdf(x,0,2);plot(x,y,x,z)

2)概率分布：P=normcdf(x,mu,sigma) 计算标准正态分布的概率P{-1

命令为：P=normcdf(1)-normcdf(-1) 结果为：P =0.6827

3)逆概率分布：x=norminv(P,mu,sigma). 即求出x ，使得P{X

1α

-

u

. 2

1α

-

u

的含义是：)1,0(~N X ,P{X<2

1α

-

u

}=2

1α

-

05.0=α时，P=0.975, =975.0u norminv(0.975)=1.96

4)随机数生成：normrnd(mu,sigma,m,n).产生m*n 阶的正态分布随机数矩阵.

4、频数直方图的描绘

1)给出数组data 的频数表的命令为： [N,X]=hist(data,k)

此命令将区间[min(data),max(data)]分为k个小区间（缺省为10），返回数组data落在每一个小区间的频数N和每一个小区间的中点X.

2)描绘数组data的频数直方图的命令为： hist(data,k)

5、参数估计

1) 设总体服从正态分布，则其点估计和区间估计可同时由以下命令获得：

[muhat,sigmahat,muci,sigmaci] = normfit(X,alpha)

此命令在显著性水平alpha下估计数据X的参数（alpha缺省时设定为0.05），返回值muhat 是X的均值的点估计值，sigmahat是标准差的点估计值, muci是均值的区间估计,sigmaci 是标准差的区间估计.

2)总体为其他分布，可以用函数[phat,pci]=mle(‘分布类型’，data,alpha)，data为数据向量，phat为未知参数的最大似然估计值，pci为置信度1-alpha的置信区间。

例如：x=poissrnd(2,1,100);[x1,x2]=mle('poiss',x,0.01)

注：

使用Matlab工具箱中具有特定分布总体的估计命令：

（a）[muhat, muci] = expfit(X,alpha)----- 在显著性水平alpha下，求指数分布的数据X的均值的点估计及其区间估计.

（b）[lambdahat, lambdaci] = poissfit(X,alpha)----- 在显著性水平alpha下，求泊松分布的数据X 的参数的点估计及其区间估计.

（c）[phat, pci] = weibfit(X,alpha)----- 在显著性水平alpha下，求Weibull分布的数据X 的参数的点估计及其区间估计.

6、假设检验

在总体服从正态分布的情况下，可用以下命令进行假设检验.

1）总体方差sigma2已知时，总体均值的检验使用 z-检验

[h,sig,ci] = ztest(x,m,sigma,alpha,tail)

检验数据 x 的关于均值的某一假设是否成立，其中sigma 为已知均方差， alpha 为显著性水平，究竟检验什么假设取决于 tail 的取值：

tail = 0，检验假设“x 的均值等于m ”

tail = 1，检验假设“x 的均值大于m ”（左侧）

tail =-1，检验假设“x 的均值小于m ”（右侧）

tail的缺省值为 0， alpha的缺省值为 0.05.

返回值 h 为一个布尔值，h=1 表示可以拒绝假设，h=0 表示接受假设，sig 为假设成立

的概率，ci 为均值的 1-alpha 置信区间.

2）总体方差sigma2未知时，总体均值的检验使用t-检验

[h,sig,ci] = ttest(x,m,alpha,tail)

检验数据 x 的关于均值的某一假设是否成立，其中alpha 为显著性水平，究竟检验什么假设取决于 tail 的取值：

tail = 0，检验假设“x 的均值等于m ”

tail = 1，检验假设“x 的均值大于m ”

tail =-1，检验假设“x 的均值小于m ”

tail的缺省值为 0， alpha的缺省值为 0.05.

返回值 h 为一个布尔值，h=1 表示可以拒绝假设，h=0 表示不可以拒绝假设，sig 为假设成立的概率，ci 为均值的 1-alpha 置信区间.

3）单正态方差的检验：

h=sctest(x,sigma2,alpha)

4）两总体均值的假设检验使用 t-检验

[h,sig,ci] = ttest2(x,y,alpha,tail)

检验数据 x ，y 的关于均值的某一假设是否成立，其中alpha 为显著性水平，究竟检验什么假设取决于 tail 的取值：

tail = 0，检验假设“x 的均值等于 y 的均值”

tail = 1，检验假设“x 的均值大于 y 的均值”

tail =-1，检验假设“x 的均值小于 y 的均值”

tail的缺省值为 0， alpha的缺省值为 0.05.

返回值 h 为一个布尔值，h=1 表示可以拒绝假设，h=0 表示不可以拒绝假设，sig 为假设成立的概率，ci 为与x与y均值差的的 1-alpha 置信区间.

注：双正态其他两个检验matlab中没有函数，所以只能按照算法自行编程实现。

5）非参数检验：总体分布的检验

主要函数youdu(p0,N),N为随机变量X的取值，即频数,p0为原假设的行向量.

注：总体分布检验没有现成的程序可以完成，需要自行编程实现，但是可以借助主要函数youdu(p0,N)来编程，这样较为方便。

程序：

alpha=0.05

p0=[1/6 1/6 1/6 1/6 1/6 1/6];

N=[21,28,19,24,16,12];

h=youdu(p0,N,alpha)

注：h=0表示接受假设，h=1拒绝原假设。

Matlab工具箱提供了两个对总体分布进行检验的命令:

（1）h = normplot(x)

此命令显示数据矩阵x的正态概率图.如果数据来自于正态分布，则图形显示出直线性形态.而其它概率分布函数显示出曲线形态.

（2）h = weibplot(x)

此命令显示数据矩阵x的Weibull概率图.如果数据来自于Weibull分布，则图形将显示出直线性形态.而其它概率分布函数将显示出曲线形态.

注：其它分布需要编程实现，matlab中没有命令。

例1某校60名学生的一次考试成绩如下:93 75 83 93 91 85 84 82 77 76 77 95 94 89 91 88 86 83 96 81 79 97 78 75 67 69 68 84 83 81 75 66 85 70 94 84 83 82 80 78 74 73 76 70 86 76 90 89 71 66 86 73 80 94 79 78 77 63 53 55

1)计算均值、标准差、极差、偏度、峰度，画出直方图；

2)检验分布的正态性；

3)若检验符合正态分布，估计正态分布的参数.

例2、在某盒中存放有白球和黑球，现作下面这样的实验：用返回抽取方式从此盒中摸球，直到摸取的是白球为止，记录下摸取的次数，重复如此的实验100次，结果见表：

找到一定方法推断盒子中的白球与黑球个数是否一样多?

程序 chengxu9

二、一元线性回归分析的Matlab实现

一元线性回归分析在matlab中用regress命令实现，用法是

b=regress(y,x)

[b,bint,r,rint,s]=regress(y,x,alpha);

输入y（因变量，列向量），X（1与自变量组成的矩阵），alpha是显著性水平。输出b系数

向量，bint 为系数置信区间，r 为残差向量，rint 为 残差置信区间（i ε为随机误差，也称残差），s 包含四个统计量（以前版本只有前三个）：第一个是样本决定系数2R ，第2个是F 值，第3个是(1,2)F n -分布大于F 的概率p ，p α<时拒绝假设，认为回归方程有效，第4个参数是2s ---剩余方差。

注：一元线性回归方程检验，可建立F 统计量，

F =

2

1?()1n

i

i y

y =-∑0

21

?()~(,2)2

n

i

i H i y

y

F p n n =---∑

若(1,2)F F n α>-，则拒绝假设，方程显著；否则接受假设，方程不显著。

rcoplot(r,rint)残差分析函数 例：

1、输入数据：

x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164]'; X=[ones(16,1) x];

Y=[88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102]'; 2、回归分析及检验：

[b,bint,r,rint,stats]=regress(Y,X)；b,bint,stats

得结果：b = -16.0730 0.7194 bint =

-33.7071 1.5612 0.6047 0.8340

stats =0.9282 180.9531 0.0000 1.7437

即0?β1?16.073,0.7194β=-=；0?β的置信区间为[-33.7017，1.5612], 1?β的置信区间为[0.6047,0.834]; r 2=0.9282, F=180.9531, p=0.0000 ，p<0.05, 可知回归模型 y=-16.073+0.7194x 成立.

例：在研究我国人均消费水平的问题中，把全国人均消费金额记作y （元），把人均国民收入

记为x （元）。我们收集到1981-1993年13年的样本数据，i=1，2，……，13，数据见表：

（1）画散点图；

（2）判断x 与y 之间是否大致成线性关系； （3）估计出回归方程；

（4）求出随机误差ε的方差2σ的估计值；

（5）给出0?β、1?β的95%的区间估计； （6）求出x 与y 的决定系数； （7）对方程进行回归分析； （8）对回归方程作残差分析；

(9)预测国民收入为2300元时人均消费金额，并给出预测区间； 解：chengxu3

三、多元线性回归

多元线性回归分析在matlab 中也主要用regress 命令实现，与一元的区别在于自变量数据 [b,bint,r,rint,s]=regress(y,x,alpha);其中X 是矩阵，而不是向量。

X =??

?

?

?

?

?

?

?np n n p p x x x x x x x x x 2

12222111211

111

利用regress 命令主要是估计参数p ββββ,,,,210 ，求出残差，以及方程显著性的检验，但是要进行具体回归分析还要用到前面所涉及的知识点计算。

下面举例说明多元分析的全过程：

例：设某商品的需求量与消费者的平均收入、商品价格的统计数据如下，建立回归模型，预测平均收入为1000、价格为6时的商品需求量.

需求量 100 75 80

70

50

65

90

100

110

60

收入

1000 600 1200 500

300

400

1300 1100 1300 300

价格 5 7 6 6 8 7 5 4 3 9

解：首先用一次函数分析：

y=[100 75 80 70 50 65 90 100 110 60];

x1=[1000 600 1200 500 300 400 1300 1100 1300 300]; x2=[5 7 6 6 8 7 5 4 3 9]; subplot(1,2,1),plot(x1,y,'g*') subplot(1,2,2),plot(x2,y,'ro') pause

X=[ones(10,1) x1' x2'];

[b,bint,r,rint,stats]=regress(y',X); b=b',bint=bint',r,rint=rint',stats 见结果:

b=111.6918 0.0143 -7.1882 bint = 56.0503 -0.0120 -13.2306 167.3334 0.0406 -1.1458

r = 9.9523 5.0477 -5.7188 -5.7109 -8.4750 -2.0929 -4.3368 1.3344 1.2867 8.7133 rint =

-4.2550 -11.3965 -17.7850 -19.9338 -22.0427 -18.1130 24.1597 21.4918 6.3474 8.5121 5.0927 13.9271 -18.5571 -14.5248 -12.6974 -2.5272 9.8836 17.1936 15.2709 19.9537 stats = 0.8944 29.6533 0.0004 52.0311 明显不太理想，所以选择二次函数

将22

011223142y x x x x βββββ=++++化为多元线性回归：

X=[ones(10,1) x1' x2' (x1.^2)' (x2.^2)']; [b,bint,r,rint,stats]=regress(y ’,X);

b,stats

结果为: b = 110.5313 0.1464 -26.5709 -0.0001 1.8475

bint = 57.2602 0.0408 -43.2247 -0.0001 0.3745

163.8024 0.2521 -9.9171 -0.0000 3.3205

r= 5.2724 -0.7162 -4.5158 -1.9390 -3.3315 3.4566 3.4843 -3.4452 -0.0976 1.8320

rint =

-2.8991 -10.7426 -11.2788 -11.3778 -12.3214 -5.9980

13.4438 9.3103 2.2472 7.4997 5.6583 12.9111

-3.5514 -13.0340 -6.3831 -3.3221

10.5200 6.1437 6.1878 6.9862

stats = 0.9702 40.6656 0.0005 20.5771

预测：y=b*[1,1000,6,1000^2,36]' 结果：y =88.4791

例：财政收入预测问题：财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据，试构造预测模型。

三、逐步回归分析

当多元线性回归方程的自变量个数不是很多时，我们有多种方法来找到‘最优’的回归方程，比如：（1）从拟合的角度考虑主要的两个准则：1.自由度调整复相关系数达到最大；2.平均残差和达到最小；（2）从极大似然估计法考虑的准则：赤池信息量AIC 达到最小；（3）从预测角度考虑的准则：统计量p C 达到最小等等方法。但是自变量,21,,,p x x x 的变量子集共有21p ，当p 很大时，以上方法的计算量是非常大的，为了解决这个问题人们提出了逐步回归的方法，帮助找到最优方程。具体过程方法前面理论部分已经详细阐述，借助于数学软件可以简化很多运算，下面用matlab 加以说明。

逐步回归的命令是： stepwise （x ，y ，[inmodel]，penter ，premove ）

x自变量数据,n*m阶矩阵；y因变量数据，n*1阶矩阵；

inmodel：矩阵的列数的指标，给出初始模型中包括的子集（缺省时设定为全部自变量）penter：引入显著性水平，默认为0.05

premove：剔除显著性水平，默认为0.1，剔除比引入要求宽松一些。

运行stepwise命令时产生一个图形窗口：列出了一个统计表，包括回归系数及其置信区间，以及模型的统计量剩余标准差（RMSE）、决定系数（R-square）、相关系数adj R（相关系数检验），F值、与F对应的概率P.

例水泥凝固时放出的热量y与水泥中4种化学成分x1、x2、x3、 x4有关，今测得一组数据如下，试用逐步回归法确定一个线性模型.

1、数据输入：

x1=[7 1 11 11 7 11 3 1 2 21 1 11 10]';

x2=[26 29 56 31 52 55 71 31 54 47 40 66 68]';

x3=[6 15 8 8 6 9 17 22 18 4 23 9 8]';

x4=[60 52 20 47 33 22 6 44 22 26 34 12 12]';

y=[78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3 109.4]';

x=[x1 x2 x3 x4];

2、逐步回归：

（1）先在初始模型中取全部自变量：stepwise(x,y)

（2）根据图表进行变量调整

（3）对变量y和x1、x2作线性回归：

得结果：b =52.5773 1.4683 0.6623,故最终模型为：y=52.5773+1.4683x1+0.6623x2 Chengxu17

注：结果并不唯一，还要结合实际情况决定。

例 调查了12名6至12岁正常儿童的体重、身高和年龄，如表所示，建立回归模型用于从身高和年龄预测儿童的体重。记儿童的体重为因变量y ，身高和年龄分别为自变量12,x x ，分别画出y 与1x 和y 与2x 的散点图，可以确定它们之间存的关系是二次关系，于是可以将

22

123142512,,,,x x x x x x x x x ===作为候选变量进行逐步回归分析，得到较优的方程。

解：chengxu5 例：

问题：一家高技术公司人事部为研究开发人员的薪金与他们的资力、管理水平、教育水平等因素之间的关系，要建立一个数学模型，以便分析人事策略的合理性，并作为新聘用人员薪金的参考。他们认为目前公司人员的薪金总体上是合理的，可以作为建模的依据，于是调查了46名软件开发人员的档案资料，其中资力一列指从事专业工作的年数，管理水平一列中1表示管理人员，0表示非管理人员，教育水平一列中1表示中学水平，2表示大学水平，3表示研究生水平。