多元回归程序MATLAB程序Word下载.doc
《多元回归程序MATLAB程序Word下载.doc》由会员分享，可在线阅读，更多相关《多元回归程序MATLAB程序Word下载.doc（17页珍藏版）》请在冰点文库上搜索。
147
149
150
153
154
155
156
157
158
159
160
162
164
腿长
88
85
91
92
93
95
96
98
97
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,'
）
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
8/30
9/30
10/30
11/30
12/30
13/30
14/30
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:
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时的商品需求量.
需求量
75
80
70
50
65
90
110
60
收入
1000
600
1200
500
300
400
1300
1100
价格
5
7
6
8
4
3
9
选择纯二次模型，即
----用户可以任意设计函数
x1=[10006001200500300400130011001300300];
x2=[5766875439];
y=[10075807050659010011060]'
X=[ones（10,1）x1'
x2'
（x1.^2）'
（x2.^2）'
];
%10指有10组数据，x1'
时方程的自变量
[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：
有n项（n个待定系数），x就有n列；
有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（:
%系数是数组，b
（1）,b
（2）,…b（n）依次代表系数1,系数2,……系数n
x2=x（:
%自变量x是一个矩阵,它的每一列分别代表一个变量,有n列就可以最多n
x3=x（:
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）;
o'
xx,yy,'
r'
）%,xx,yy,'
是画曲线，相当于拟合
结果：
a=87.5244
k1=0.0269
k2=-63.4591
m=0.1083
图形：
例3
出钢时所用的盛钢水的钢包，由于钢水对耐火材料的侵蚀，容积不断增大.我们希望知道使用次数与增大的容积之间的关系.对一钢包作试验，测得的数据列于下表：
使用次数
增大容积
2
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.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];
*-'
beta0=[82]'
[beta,r,J]=nlinfit（x'
y'
zhang1'
beta0）
3、求回归系数：
beta=
11.6037
-1.0641
即得回归模型为：
ro'
）
holdon
xx=2:
0.05:
yy=zhang1（beta,xx）;
%-通过调用用户自编的函数
plot（xx,yy,'
）%拟合成线
或者
[YY,delta]=nlpredci（'
x'
beta,r,J）
k+'
x,YY,'
yy=beta
（1）*exp（beta
（2）./xx）;
plot（xx,yy,'
例4
财政收入预测问题：
财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。
下表列出了1952-1981年的原始数据，试构造预测模型。
财政收入预测问题：
年份
国民收入（亿元）
工业总产值（亿元）
农业总产值（亿元）
总人口（万人）
就业人口（万人）
固定资产投资（亿元）
财政收入（亿元）
1952
598
349
461
57482
20729
44
184
1953
586
455
475
58796
21364
89
216
1954
707
520
491
60266
21832
248
1955
737
558
529
61465
22328
254
1956
825
715
556
62828
23018
268
1957
837
798
575
64653
23711
139
286
1958
1028
1235
65994
26600
256
357
1959
1114
1681
509
67207
26173
338
444
1960
1079
1870
66207
25880
380
506
1961
757
1156
434
65859
25590
138
271
1962
677
964
67295
25110
66
230
1963
779
1046
514
69172
26640
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
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
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
2475
5590
1058
96259
39856
550
922
1979
2702
6065
1150
97542
40581
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）%分别代表每个系数，而第二个参数代表所有的自变量，%是一个矩阵，它的每一列分别代表一个自变量。
升级会员 