Matlab统计工具箱的应用.docx

Matlab统计工具箱的应用.docx

《Matlab统计工具箱的应用.docx》由会员分享，可在线阅读，更多相关《Matlab统计工具箱的应用.docx（19页珍藏版）》请在冰豆网上搜索。

Matlab统计工具箱的应用

Matlab统计工具箱的应用

一、统计的基本命令

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

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

年份

78

79

80

81

82

82

84

85

86

87

职工工资总额（亿元）

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

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

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

命令格式：

x=a:

h:

bt=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中:

savedatatxy

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

loaddata

方法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中：

savedata1data

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

loaddata1

调用数据文件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）loaddata.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

分布：

chi2t分布：

tF分布：

F

Matlab工具箱对每一种分布都提供五类函数，其命令字符为：

概率密度：

pdf概率分布：

cdf逆概率分布：

inv

均值与方差：

stat随机数生成：

rnd

（当需要一种分布的某一类函数时，将以上所列的分布命令字符与函数命令字符接起来，并输入自变量（可以是标量、数组或矩阵）和参数即可.

如对均值为mu、标准差为sigma的正态分布，举例如下：

例画出正态分布

和

的概率密度函数图形.

在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

例取

，求

.

的含义是：

P{X<

}=

时，P=0.975,

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）来编程，这样较为方便。

例如：

投色子120次统计如下频数

点数

123456

频数

212819241612

问色子是否均匀？

（

）

程序：

alpha=0.05

p0=[1/61/61/61/61/61/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名学生的一次考试成绩如下:

937583939185848277767795948991888683968179977875676968848381756685709484838280787473767086769089716686738094797877635355

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

2）检验分布的正态性；

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

例2、在某盒中存放有白球和黑球，现作下面这样的实验：

用返回抽取方式从此盒中摸球，直到摸取的是白球为止，记录下摸取的次数，重复如此的实验100次，结果见表：

摸取次数

1234

频数

43311565

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

程序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为残差置信区间（

为随机误差，也称残差），s包含四个统计量（以前版本只有前三个）：

第一个是样本决定系数

，第2个是

值，第3个是

分布大于

的概率p，

时拒绝假设，认为回归方程有效，第4个参数是

---剩余方差。

注：

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

若

，则拒绝假设，方程显著；否则接受假设，方程不显著。

rcoplot（r,rint）残差分析函数

例：

1、输入数据：

x=[143145146147149150153154155156157158159160162164]';

X=[ones（16,1）x];

Y=[8885889192939395969897969899100102]';

2、回归分析及检验：

[b,bint,r,rint,stats]=regress（Y,X）；b,bint,stats

得结果：

b=

-16.0730

0.7194

bint=

-33.70711.5612

0.60470.8340

stats=0.9282180.95310.00001.7437

即

；

的置信区间为[-33.7017，1.5612],

的置信区间为[0.6047,0.834];r2=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，数据见表：

年份

人均国民收入x（元）

人均消费金额y（元）

年份

人均国民收入x（元）

人均消费金额y（元）

1981

1982

1983

1984

1985

1986

1987

393.8

419.14

460.86

544.11

668.29

737.73

859.97

249

267

289

329

406

451

513

1988

1989

1990

1991

1992

1993

1068.8

1169.2

1250.7

1429.5

1725.9

2099.5

643

699

713

803

947

1148

（1）画散点图；

（2）判断x与y之间是否大致成线性关系；

（3）估计出回归方程；

（4）求出随机误差

的方差

的估计值；

（5）给出

、

的95%的区间估计；

（6）求出x与y的决定系数；

（7）对方程进行回归分析；

（8）对回归方程作残差分析；

（9）预测国民收入为2300元时人均消费金额，并给出预测区间；

解：

chengxu3

三、多元线性回归

多元线性回归分析在matlab中也主要用regress命令实现，与一元的区别在于自变量数据

[b,bint,r,rint,s]=regress（y,x,alpha）;其中X是矩阵，而不是向量。

X=

利用regress命令主要是估计参数

，求出残差，以及方程显著性的检验，但是要进行具体回归分析还要用到前面所涉及的知识点计算。

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

例：

设某商品的需求量与消费者的平均收入、商品价格的统计数据如下，建立回归模型，预测平均收入为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=[10075807050659010011060];

x1=[10006001200500300400130011001300300];

x2=[5766875439];

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.69180.0143-7.1882

bint=56.0503-0.0120-13.2306

167.33340.0406-1.1458

r=9.95235.0477-5.7188-5.7109-8.4750-2.0929-4.33681.33441.28678.7133

rint=

-4.2550-11.3965-17.7850-19.9338-22.0427-18.1130

24.159721.49186.34748.51215.092713.9271

-18.5571-14.5248-12.6974-2.5272

9.883617.193615.270919.9537

stats=0.894429.65330.000452.0311

明显不太理想，所以选择二次函数

将

化为多元线性回归：

X=[ones（10,1）x1'x2'（x1.^2）'（x2.^2）'];

[b,bint,r,rint,stats]=regress（y’,X）;

b,stats

结果为:

b=110.53130.1464-26.5709-0.00011.8475

bint=57.26020.0408-43.2247-0.00010.3745

163.80240.2521-9.9171-0.00003.3205

r=5.2724-0.7162-4.5158-1.9390-3.33153.45663.4843-3.4452-0.09761.8320

rint=

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

13.44389.31032.24727.49975.658312.9111

-3.5514-13.0340-6.3831-3.3221

10.52006.14376.18786.9862

stats=0.970240.66560.000520.5771

预测：

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

y=88.4791

例：

财政收入预测问题：

财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。

下表列出了1952-1981年的原始数据，试构造预测模型。

年份

国民收入（亿元）

工业总产值（亿元）

农业总产值（亿元）

总人口（万人）

就业人口（万人）

固定资产投资（亿元）

财政收入（亿元）

1952

598

349

461

57482

20729

44

184

1953

586

455

475

58796

21364

89

216

1954

707

520

491

60266

21832

97

248

1955

737

558

529

61465

22328

98

254

1956

825

715

556

62828

23018

150

268

1957

837

798

575

64653

23711

139

286

1958

1028

1235

598

65994

26600

256

357

1959

1114

1681

509

67207

26173

338

444

1960

1079

1870

444

66207

25880

380

506

1961

757

1156

434

65859

25590

138

271

1962

677

964

461

67295

25110

66

230

1963

779

1046

514

69172

26640

85

266

1964

943

1250

584

70499

27736

129

323

1965

1152

1581

632

72538

28670

175

393

1966

1322

1911

687

74542

29805

212

466

1967

1249

1647

697

76368

30814

156

352

1968

1187

1565

680

78534

31915

127

303

1969

1372

2101

688

80671

33225

207

447

1970

1638

2747

767

82992

34432

312

564

1971

1780

3156

790

85229

35620

355

638

1972

1833

3365

789

87177

35854

354

658

1973

1978

3684

855

89211

36652

374

691

1974

1993

3696

891

90859

37369

393

655

1975

2121

4254

932

92421

38168

462

692

1976

2052

4309

955

93717

38834

443

657

1977

2189

4925

971

94974

39377

454

723

1978

2475

5590

1058

96259

39856

550

922

1979

2702

6065

1150

97542

40581

564

890

1980

2791

6592

1194

98705

41896

568

826

1981

2927

6862

1273

100072

73280

496

810

三、逐步回归分析

当多元线性回归方程的自变量个数不是很多时，我们有多种方法来找到‘最优’的回归方程，比如：

（1）从拟合的角度考虑主要的两个准则：

1.自由度调整复相关系数达到最大；2.平均残差和达到最小；

（2）从极大似然估计法考虑的准则：

赤池信息量AIC达到最小；（3）从预测角度考虑的准则：

统计量

达到最小等等方法。

但是自变量

的变量子集共有

，当p很大时，以上方法的计算量是非常大的，为了解决这个问题人们提出了逐步回归的方法，帮助找到最优方程。

具体过程方法前面理论部分已经详细阐述，借助于数学软件可以简化很多运算，下面用matlab加以说明。

逐步回归的命令是：

stepwise（x，y，[inmodel]，penter，premove）

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

inmodel：

矩阵的列数的指标，给出初始模型中包括的子集（缺省时设定为全部自变量）

penter：

引入显著性水平，默认为0.05

premove：

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

运行stepwi