多元回归程序MATLAB程序.doc

上传人:wj 文档编号:2501609 上传时间:2023-05-03 格式:DOC 页数:17 大小:208.50KB
下载 相关 举报
多元回归程序MATLAB程序.doc_第1页
第1页 / 共17页
多元回归程序MATLAB程序.doc_第2页
第2页 / 共17页
多元回归程序MATLAB程序.doc_第3页
第3页 / 共17页
多元回归程序MATLAB程序.doc_第4页
第4页 / 共17页
多元回归程序MATLAB程序.doc_第5页
第5页 / 共17页
多元回归程序MATLAB程序.doc_第6页
第6页 / 共17页
多元回归程序MATLAB程序.doc_第7页
第7页 / 共17页
多元回归程序MATLAB程序.doc_第8页
第8页 / 共17页
多元回归程序MATLAB程序.doc_第9页
第9页 / 共17页
多元回归程序MATLAB程序.doc_第10页
第10页 / 共17页
多元回归程序MATLAB程序.doc_第11页
第11页 / 共17页
多元回归程序MATLAB程序.doc_第12页
第12页 / 共17页
多元回归程序MATLAB程序.doc_第13页
第13页 / 共17页
多元回归程序MATLAB程序.doc_第14页
第14页 / 共17页
多元回归程序MATLAB程序.doc_第15页
第15页 / 共17页
多元回归程序MATLAB程序.doc_第16页
第16页 / 共17页
多元回归程序MATLAB程序.doc_第17页
第17页 / 共17页
亲,该文档总共17页,全部预览完了,如果喜欢就下载吧!
下载资源
资源描述

多元回归程序MATLAB程序.doc

《多元回归程序MATLAB程序.doc》由会员分享,可在线阅读,更多相关《多元回归程序MATLAB程序.doc(17页珍藏版)》请在冰点文库上搜索。

多元回归程序MATLAB程序.doc

matlab回归(拟合)总结

前言

1、学三条命令

polyfit(x,y,n)---拟合成一元幂函数(一元多次)

regress(y,x)----可以多元,

nlinfit(x,y,’fun’,beta0)(可用于任何类型的函数,任意多元函数,应用范围最广,最万能的)

2、同一个问题,这三条命令都可以使用,但结果肯定是不同的,因为拟合的近似结果,没有唯一的标准的答案。

相当于咨询多个专家。

3、回归的操作步骤:

根据图形(实际点),选配一条恰当的函数形式(类型)---需要数学理论与基础和经验。

(并写出该函数表达式的一般形式,含待定系数)------选用某条回归命令求出所有的待定系数。

所以可以说,回归就是求待定系数的过程(需确定函数的形式)

一、多元回归分析

对于多元线性回归模型(其实可以是非线性,它通用性极高):

设变量的n组观测值为

记,,则的估计值为排列方式与线性代数中的线性方程组相同(),拟合成多元函数---regress

使用格式:

左边用b=[b,bint,r,rint,stats]右边用=regress(y,x)或regress(y,x,alpha)

---命令中是先y后x,

---须构造好矩阵x(x中的每列与目标函数的一项对应)

---并且x要在最前面额外添加全1列/对应于常数项

---y必须是列向量

---结果是从常数项开始---与polyfit的不同。

其中:

b为回归系数,

的估计值(第一个为常数项),

bint为回归系数的区间估计,

r:

残差,

rint:

残差的置信区间,

stats:

用于检验回归模型的统计量,有四个数值:

相关系数r2、F值、与F对应的概率p和残差的方差(前两个越大越好,后两个越小越好),

alpha:

显著性水平(缺省时为0.05,即置信水平为95%),(alpha不影响b,只影响bint(区间估计)。

它越小,即置信度越高,则bint范围越大。

显著水平越高,则区间就越小)(返回五个结果)---

如有n个自变量-有误(n个待定系数),则b中就有n+1个系数(含常数项,---第一项为常数项)(b---b的范围/置信区间---残差r---r的置信区间rint-----点估计----区间估计

如果的置信区间(bint的第行)不包含0,则在显著水平为时拒绝的假设,认为变量是显著的.*******(而rint残差的区间应包含0则更好)。

b,y等均为列向量,x为矩阵(表示了一组实际的数据)必须在x第一列添加一个全1列。

----对应于常数项。

相关系数r2越接近1,说明回归方程越显著;(r2越大越接近1越好)F越大,说明回归方程越显著;(F越大越好)与F对应的概率p越小越好,一定要P

乘余(残差)标准差(RMSE)越小越好(此处是残差的方差,还没有开方)(前两个越大越好,后两个越小越好)

重点:

regress(y,x)重点与难点是如何加工处理矩阵x。

y是函数值,一定是只有一列。

也即目标函数的形式是由矩阵X来确定

如s=a+b*x1+c*x2+d*x3+e*x1^2+f*x2*x3+g*x1^2,

一定有一个常数项,且必须放在最前面(即x的第一列为全1列)

X中的每一列对应于目标函数中的一项(目标函数有多少项则x中就有多少列)

X=[ones,x1,x2,x3,x1.^2,x2.*x3,x1.ˆ2](剔除待定系数的形式)

regress:

y/x顺序,矩阵X需要加工处理

nlinfit:

x/y顺序,X/Y就是原始的数据,不要做任何的加工。

(即regress靠矩阵X来确定目标函数的类型形式(所以X很复杂,要作很多处理)而nlinfit是靠程序来确定目标函数的类型形式(所以X就是原始数据,不要做任何处理)

例1

测16名成年女子的身高与腿长所得数据如下:

身高

143

145

146

147

149

150

153

154

155

156

157

158

159

160

162

164

腿长

88

85

88

91

92

93

93

95

96

98

97

96

98

99

100

102

配成y=a+b*x形式

x=[143145146147149150153154155156157158159160162164]';

y=[8885889192939395969897969899100102]';

plot(x,y,'r+')

z=x;

x=[ones(16,1),x];----常数项

[b,bint,r,rint,stats]=regress(y,x);---处结果与polyfit(x,y,1)相同

b,bint,stats

得结果:

b=bint=

-16.0730-33.70711.5612------每一行为一个区间

0.71940.60470.8340

stats=0.9282180.95310.0000

即;的置信区间为[-33.7017,1.5612],的置信区间为[0.6047,0.834];r2=0.9282,F=180.9531,p=0.0。

p<0.05,可知回归模型y=-16.073+0.7194x成立.

[b,bint,r,rint,stats]=regress(Y,X,0.05);-----结果相同

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

polyfit(x,y,1)-----当为一元时(也只有一组数),则结果与regress是相同的,只是

命令中x,y要交换顺序,结果的系数排列顺序完全相反,x中不需要全1列。

ans=0.7194-16.0730--此题也可用polyfit求解,杀鸡用牛刀,脖子被切断。

3、残差分析,作残差图:

rcoplot(r,rint)

从残差图可以看出,除第二个数据外,其余数据的残差离零点均较近,且残差的置信区间均包含零点,这说明回归模型y=-16.073+0.7194x能较好的符合原始数据,而第二个数据可视为异常点(而剔除)

4、预测及作图:

plot(x,y,'r+')

holdon

a=140:

165;

b=b

(1)+b

(2)*a;

plot(a,b,'g')

例2

观测物体降落的距离s与时间t的关系,得到数据如下表,求s关于t的回归方程

t(s)

1/30

2/30

3/30

4/30

5/30

6/30

7/30

s(cm)

11.86

15.67

20.60

26.69

33.71

41.93

51.13

t(s)

8/30

9/30

10/30

11/30

12/30

13/30

14/30

s(cm)

61.49

72.90

85.44

99.08

113.77

129.54

146.48

法一:

直接作二次多项式回归

t=1/30:

1/30:

14/30;

s=[11.8615.6720.6026.6933.7141.9351.1361.4972.9085.4499.08113.77129.54146.48];

[p,S]=polyfit(t,s,2)

p=489.294665.88969.1329

得回归模型为:

方法二----化为多元线性回归:

t=1/30:

1/30:

14/30;

s=[11.8615.6720.6026.6933.7141.9351.1361.4972.9085.4499.08113.77129.54146.48];

T=[ones(14,1),t',(t.^2)']%?

是否可行?

等验证...----因为有三个待定系数,所以有三列,始于常数项

[b,bint,r,rint,stats]=regress(s',T);

b,stats

b=9.1329

65.8896

489.2946

stats=1.0e+007*

0.00001.037800.0000

得回归模型为:

%结果与方法1相同

polyfit------一元多次

regress----多元一次---其实通过技巧也可以多元多次

regress最通用的,万能的,表面上是多元一次,其实可以变为多元多次且任意函数,如x有n列(不含全1列),则表达式中就有n+1列(第一个为常数项,其他每项与x的列序相对应)。

例3

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

选择纯二次模型,即

----用户可以任意设计函数

x1=[10006001200500300400130011001300300];

x2=[5766875439];

y=[10075807050659010011060]';

X=[ones(10,1)x1'x2'(x1.^2)'(x2.^2)'];%10指有10组数据,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

stats=0.970240.66560.000520.5771

故回归模型为:

剩余标准差为4.5362,说明此回归模型的显著性较好.

三、非线性回归(拟合)

使用格式:

beta=nlinfit(x,y,‘程序名’,beta0)

[beta,r,J]=nlinfit(X,y,fun,beta0)

X给定的自变量数据,

Y给定的因变量数据,

fun要拟合的函数模型(句柄函数或者内联函数形式),

beta0函数模型中待定系数估计初值(即程序的初始实参)

beta返回拟合后的待定系数

其中beta为估计出的回归系数;r为残差;J为Jacobian矩阵

输入数据x、y分别为n*m矩阵和n维列向量,对一元非线性回归,x为n维列向量。

可以拟合成任意函数。

最通用的,万能的命令

x,y顺序,x不需要任何加工,直接用原始数据。

---所编的程序一定是两个形参(待定系数/向量,自变量/矩阵:

每一列为一个自变量)

结果要看残差的大小和是否有警告信息,如有警告则换一个b0初始向量再重新计算。

本程序中也可能要用.*./.^如结果中有警告信息,则必须多次换初值来试算.

难点是编程序与初值

存在的问题:

不同的beta0,则会产生不同的结果,如何给待定系数的初值以及如何分析结果的好坏,如出现警告信息,则换一个待定系数试一试。

因为拟合本来就是近似的,可能有多个结果。

1:

重点(难点)是预先编程序(即确定目标函数的形式,而regress的目标函数由x矩阵来确定,其重难点为构造矩阵a)

2:

x/y顺序—列向量----x/y是原始数据,不要做任何修改

3:

编程:

一定两个形参(beta,x)a=beta

(1);b=beta

(2);c=beta(3);…x1=x(:

1);x2=x(:

2);x3=x(:

3);即每一列为一个自变量

4:

regress/nlinfit都是列向量

5:

regress:

有n项(n个待定系数),x就有n列;nlinfit:

有m个变量则x就有m列

例1

已知数据:

x1=[0.5,0.4,0.3,0.2,0.1]; x2=[0.3,0.5,0.2,0.4,0.6];

x3=[1.8,1.4,1.0,1.4,1.8];y=[0.785,0.703,0.583,0.571,0.126]’;且y与x1,x2,x3关系为多元非线性关系(只与x2,x3相关)为:

 y=a+b*x2+c*x3+d*(x2.^2)+e*(x3.^2)—此函数是由用户根据图形的形状等所配的曲线,即自己选定函数类型求非线性回归系数a,b,c,d,e。

(1)对回归模型建立M文件model.m如下:

functionyy=myfun(beta,x)%一定是两个参数:

系数和自变量---一个向量/一个矩阵

a=beta

(1)

b=beta

(2)

c=beta(3)

d=beta(4)

e=beta(5)

x1=x(:

1);%系数是数组,b

(1),b

(2),…b(n)依次代表系数1,系数2,……系数n

x2=x(:

2);%自变量x是一个矩阵,它的每一列分别代表一个变量,有n列就可以最多n

x3=x(:

3);

yy=beta

(1)+beta

(2)*x2+beta(3)*x3+beta(4)*(x2.^2)+beta(5)*(x3.^2);

(b(i)与待定系数的顺序关系可以任意排列,并不是一定常数项在最前,只是结果与自己指定的相对应)(x一定是一列对应一个变量,不能x1=x

(1),x2=x

(2),x3=x(3)……)

(2)主程序如下:

x=[0.5,0.4,0.3,0.2,0.1;0.3,0.5,0.2,0.4,0.6;1.8,1.4,1.0,1.4,1.8]';%每一列为一个变量,如果是倒入数据矩阵,只能把x的数据倒进去,不能把全部数据都倒进去,然后选某几列

y=[0.785,0.703,0.583,0.571,0.126]';

beta0=[1,1,1,1,1,1]';%有多少个待定系数,就给多少个初始值。

[beta,r,j]=nlinfit(x,y,@myfun,beta0)

beta=-0.44205.51110.3837-8.1734-0.1340

例2

混凝土的抗压强度随养护时间的延长而增加,现将一批混凝土作成12个试块,记录了养护日期(日)及抗压强度y(kg/cm2)的数据:

养护时间:

x=[234579121417212856]抗压强度:

y=[35+r42+r47+r53+r59+r65+r68+r73+r76+r82+r86+r99+r]建立非线性回归模型,对得到的模型和系数进行检验。

注明:

此题中的+r代表加上一个[-0.5,0.5]之间的随机数模型为:

y=a+k1*exp(m*x)+k2*exp(-m*x);------有四个待定系数

Matlab程序:

x=[234579121417212856];

r=rand(1,12)-0.5;

y1=[354247535965687376828699];

y=y1+r;

myfunc=inline('beta

(1)+beta

(2)*exp(beta(4)*x)+beta(3)*exp(-beta(4)*x)','beta','x');

beta=nlinfit(x,y,myfunc,[0.50.50.50.5]);%初值为0.2也可以,如为1则不行,则试着换系数初值----此处为一元,x’,y’行/列向量都可以

a=beta

(1),k1=beta

(2),k2=beta(3),m=beta(4)  %testthemodel

xx=min(x):

max(x); %2:

56

yy=a+k1*exp(m*xx)+k2*exp(-m*xx);

plot(x,y,'o',xx,yy,'r')%,xx,yy,'r'是画曲线,相当于拟合

结果:

a=87.5244

k1=0.0269

k2=-63.4591

m=0.1083

图形:

例3

出钢时所用的盛钢水的钢包,由于钢水对耐火材料的侵蚀,容积不断增大.我们希望知道使用次数与增大的容积之间的关系.对一钢包作试验,测得的数据列于下表:

使用次数

增大容积

使用次数

增大容积

2

3

4

5

6

7

8

9

6.42

8.20

9.58

9.50

9.70

10.00

9.93

9.99

10

11

12

13

14

15

16

10.49

10.59

10.60

10.80

10.60

10.90

10.76

对将要拟合的非线性模型y=aeb/x,(如再加y=c*sin(x)+aeb/x)

建立m-文件volum.m如下:

functionyhat=volum(beta,x)

yhat=beta

(1)*exp(beta

(2)./x);或

functionf=volum(beta,x)

a=beta

(1);

b=beta

(2);

f=a*exp(b./x);

2、输入数据:

主程序:

x=2:

16;

y=[6.428.209.589.59.7109.939.9910.4910.5910.6010.8010.6010.9010.76];

plot(x,y,'*-')

beta0=[82]';

[beta,r,J]=nlinfit(x',y','zhang1',beta0)

3、求回归系数:

beta=

11.6037

-1.0641

即得回归模型为:

4、预测及作图:

plot(x,y,'ro')

holdon

xx=2:

0.05:

16;

yy=zhang1(beta,xx);%-通过调用用户自编的函数

plot(xx,yy,'g')%拟合成线

或者

[YY,delta]=nlpredci('zhang1',x',beta,r,J)

plot(x,y,'k+',x,YY,'r')

plot(x,y,'ro')

holdon

xx=2:

0.05:

16;

yy=beta

(1)*exp(beta

(2)./xx);

plot(xx,yy,'g')

例4

财政收入预测问题:

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

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

财政收入预测问题:

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

下表列出了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

解设国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资分别为x1、x2、x3、x4、x5、x6,财政收入为y,设变量之间的关系为:

y=ax1+bx2+cx3+dx4+ex5+fx6

使用非线性回归方法求解。

1. 对回归模型建立M文件model.m如下:

functionyy=model(beta0,X)%一定是两个参数,第一个为系数数组,b

(1),b

(2),…b(n)%分别代表每个系数,而第二个参数代表所有的自变量,%是一个矩阵,它的每一列分别代表一个自变量。

展开阅读全文
相关资源
猜你喜欢
相关搜索
资源标签

当前位置:首页 > 党团工作 > 入党转正申请

copyright@ 2008-2023 冰点文库 网站版权所有

经营许可证编号:鄂ICP备19020893号-2