1:学三条命令:polyfit(x,y,n)---拟合成一元幂函数(一元多次)
regress(y,x)----可以多元, nlinfit(x,y,’fun’,beta0) (可用于任何类型的函数,任意多元函数,应用范围最主,最万能的)
2:同一个问题,可能这三条命令都可以使用,但结果肯定是不同的,因为拟合的近似结果,没有唯一的标准的答案。相当于咨询多个专家。
3:回归的操作步骤:
(1)根据图形(实际点),选配一条恰当的函数形式(类型)---需要数学理论与基础和经验。(并写出该函数表达式的一般形式,含待定系数)
(2)选用某条回归命令求出所有的待定系数
所以可以说,回归就是求待定系数的过程(需确定函数的形式)
配曲线的一般方法是:
(一)先对两个变量x和y 作n次试验观察得画出散点图,
散点图
(二)根据散点图确定须配曲线的类型.
通常选择的六类曲线如下:
(1)双曲线
(2)幂函数曲线y=a, 其中x>0,a>0
(3)指数曲线y=a其中参数a>0.
(4)倒指数曲线y=a其中a>0,
(5)对数曲线y=a+blogx,x>0
(6)S型曲线
(三)然后由n对试验数据确定每一类曲线的未知参数a和b.
一、一元多次拟合polyfit(x,y,n)
一元回归polyfit
多元回归regress---nlinfit(非线性)
二、多元回归分析
(其实可以是非线性,它通用性极高)
对于多元线性回归模型:
设变量的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列。----对应于常数项-------而nlinfit不能额外添加全1列。
结果的系数就是与此矩阵相对应的(常数项,x1,x2,……xn)。
(结果与参数个数:1/5=2/3-----y,x顺序---x要额外添加全1列)
而nlinfit:1/3=4------x,y顺序---x不能额外添加全1列,---需编程序,用于模仿需拟合的函数的任意形式,一定两个参数,一为系数数组,二为自变量矩阵(每列为一个自变量)
有n个变量---不准确,x中就有n列,再添加一个全1列(相当于常数项),就变为n+1列,则结果中就有n+1个系数。
x需要经过加工,如添加全1列,可能还要添加其他需要的变换数据。
相关系数r2越接近1,说明回归方程越显著;(r2越大越接近1越好)
F越大,说明回归方程越显著;(F越大越好)
身高 | 143 | 145 | 146 | 147 | 149 | 150 | 153 | 154 | 155 | 156 | 157 | 158 | 159 | 160 | 162 | 1 |
腿长 | 88 | 85 | 88 | 91 | 92 | 93 | 93 | 95 | 96 | 98 | 97 | 96 | 98 | 99 | 100 | 102 |
配成y=a+b*x形式
>> x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 1]';
>> y=[88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102]';
>> 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.7071 1.5612------每一行为一个区间
0.7194 0.6047 0.8340
stats =
0.9282 180.9531 0.0000
即;的置信区间为[-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 成立.
>> [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 =----此题也可用polyfit求解,杀鸡用牛刀,脖子被切断。---用牛刀杀羊
0.7194 -16.0730
3、残差分析,作残差图:
>>rcoplot(r,rint)
从残差图可以看出,除第二个数据外,其余数据的残差离零点均较近,且残差的置信区间均包含零点,这说明回归模型 y=-16.073+0.7194x能较好的符合原始数据,而第二个数据可视为异常点(而剔除)
4、预测及作图:
>> plot(x,y,'r+') >> hold on
>> 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.86 15.67 20.60 26.69 33.71 41.93 51.13 61.49 72.90 85.44 99.08 113.77 129.54 146.48];
>> [p,S]=polyfit(t,s,2)
p =
4.2946 65.86 9.1329
S =
R: [3x3 double]
df: 11
normr: 0.1157
得回归模型为 :
方法二----化为多元线性回归:
t=1/30:1/30:14/30;
s=[11.86 15.67 20.60 26.69 33.71 41.93 51.13 61.49 72.90 85.44 99.08 113.77 129.54 146.48];
>> T=[ones(14,1), t', (t.^2)'] %是否可行等验证...
T = 1.0000 0.0333 0.0011-------因为有三个待定系数,所以有三列,始于常数项
1.0000 0.0667 0.0044
1.0000 0.1000 0.0100
1.0000 0.1333 0.0178
1.0000 0.1667 0.0278
1.0000 0.2000 0.0400
1.0000 0.2333 0.0544
1.0000 0.2667 0.0711
1.0000 0.3000 0.0900
1.0000 0.3333 0.1111
1.0000 0.3667 0.1344
1.0000 0.4000 0.1600
1.0000 0.4333 0.1878
1.0000 0.4667 0.2178;
>> [b,bint,r,rint,stats]=regress(s',T);
>> b,stats
b = 9.1329
65.86
4.2946
stats =1.0e+007 *
0.0000 1.0378 0 0.0000
得回归模型为 :
%结果与方法1相同
>> T=[ones(14,1),t, (t.^2)'] %是否可行等验证...
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=[1000 600 1200 500 300 400 1300 1100 1300 300];
>> x2=[5 7 6 6 8 7 5 4 3 9];
>> y=[100 75 80 70 50 65 90 100 110 60]';
X=[ones(10,1) x1' x2' (x1.^2)' (x2.^2)']; %注意技巧性
[b,bint,r,rint,stats]=regress(y,X); %这是万能方法需进一步验证
>> b,stats
b =
110.5313
0.14
-26.5709
-0.0001
1.8475
stats = 0.9702 40.6656 0.0005 20.5771
故回归模型为:
剩余标准差为4.5362, 说明此回归模型的显著性较好.
--------(此题还可以用 rstool(X,Y)命令求解,详见回归问题详解)
>> X=[ones(10,1) x1' x2' (x1.^2)' (x2.^2)',sin(x1.*x2)',(x1.*exp(x2))'];
>> [b,bint,r,rint,stats]=regress(y,X);
>> b,stats
其他例题可以见下面的例题或练习……
(个人2011年认为,regress只能用于函数中的每一项只能有一个待定系数的情况,不能用于aebx等的情况)
2011的理解:
1/5----2/3
regress(y,x)
----re是y/x逆置的
---y是列向量
---须确定目标函数的形式
---x须构造(通过构造来反映目标函数)
---x中的每一列与目标函数的一项对应(剔除待定系数)
----首项为常数项(x的第一列为全1)
----有函数有n 项(待定系数),则x就有n列
----regress只能解决每项只有一个待定系数的情况且必须有常数项的情况(且每项只有一个待定系数,即项数与待定系数数目相同)
***其重(难、关键)点:列向量、构造矩阵(X):目标函数中的每项与X中的一列对应。(由X来确定目标函数的类型/形式)
三、非线性回归(拟合)
使用格式:
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维列向量。
’model’为是事先用m-文件定义的非线性函数;beta0为回归系数的初值
可以拟合成任意函数。最通用的,万能的命令
x,y顺序,x不需要任何加工,直接用原始数据。(也不需要全1列)---所编的程序一定是两个形参(待定系数/向量,自变量/矩阵:每一列为一个自变量)
结果要看残差的大小和是否有警告信息,如有警告则换一个b0初始向量再重新计算。
本命令的重点与难点是“编好程序”(用于确定拟拟合成的函数形式)
本程序中也可能要用.* ./ .^
如结果中有警告信息,则必须多次换初值来试算
难点是编程序与初值
自己总结:nlinfit
多元任意函数,(自己任意设计函数,再求待定系数)
13/4
顺序([b,r,j]=nlinfit(x,y,’…’, b0)
y为列向量;x为矩阵,无需加全1列,x,y就是原始的数据点,(x/y正顺序,所以x不要加全1列)
需预先编程(两个参数,系数向量,各变量的矩阵/每列为一个变量)
存在的问题:不同的beta0,则会产生不同的结果,如何给待定系数的初值以及如何分析结果的好坏
如出现警告信息,则换一个待定系数试一试。
因为拟合本来就是近似的,可能有多个结果。
2011年的理解:
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如下:
function yy=myfun(beta,x) %一定是两个参数:系数和自变量---一个向量/一个矩阵
a=beta(1)
b=beta(2)
c=beta(3)
………………………………
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]';-----每一列为一个变量
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.4420 5.5111 0.3837 -8.1734 -0.1340
此题也可用regress来求解,但结果是不一样的
>> 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]';
>> n=length(x1);
>> x=[ones(n,1),x2',x3',(x2.^2)',(x3.^2)'];
>> [b,bint,r,rint,stats]=regress(y,x);
>> b,stats
b =
-3.3844
-1.8450
6.5137
0
-2.1773
stats =
0.7859 1.2232 0.5674 0.0557
2011年题目改为:y=a+b*x1+c*x2+d*(x3.^2)+e*(x1.^2)+f*sin(x2)
求非线性回归系数a , b , c , d , e,f
function f=fxxnh(beta,x) %所编的程序一定是两个形参,第一个为待定系数向量,第二个为自变量矩阵
a=beta(1);
b=beta(2);
c=beta(3);
d=beta(4);
e=beta(5);
f=beta(6);%系数向量中的一个元素代表一个待定系数
x1=x(:,1);%自变量矩阵每一列代表一个自变量
x2=x(:,2);
x3=x(:,3);
f=a+b.*x1+c.*x2+d.*(x3.^2)+e.*(x1.^2)+f.*sin(x2);
但计算出现了问题
例题2:混凝土的抗压强度随养护时间的延长而增加,现将一批混凝土作成12个试块,记录了养护日期(日)及抗压强度y(kg/cm2)的数据:
养护时间:x =[2 3 4 5 7 9 12 14 17 21 28 56 ]
抗压强度:y =[35+r 42+r 47+r 53+r 59+r 65+r 68+r 73+r 76+r 82+r 86+r 99+r ]
建立非线性回归模型,对得到的模型和系数进行检验。
注明:此题中的+r代表加上一个[-0.5,0.5]之间的随机数
模型为:y=a+k1*exp(m*x)+k2*exp(-m*x);
---------有四个待定系数
Matlab程序:
x=[2 3 4 5 7 9 12 14 17 21 28 56];
r=rand(1,12)-0.5;
y1=[35 42 47 53 59 65 68 73 76 82 86 99];
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.5 0.5 0.5 0.5]); ---初值为0.2也可以,如为1则不行,则试着换系数初值----此处为一元,x’,y’行/列向量都可以
a=beta(1),k1=beta(2),k2=beta(3),m=beta(4)
%test the model
xx=min(x):max(x); -----2:56
yy=a+k1*exp(m*xx)+k2*exp(-m*xx);
plot(x,y,'o',xx,yy,'r')
结果:
a = 87.5244
k1 = 0.0269
k2 = -63.4591
m = 0.1083
图形:
本人认为:此题不能用regress求解,因为有些式子中含有两个待定系数
例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= c*sin(x)+aeb/x)
建立m-文件volum.m如下:
function yhat=volum(beta,x) %此函数一定是两个参数,第一个为系数数组,第二个为自变量矩阵。
yhat=beta(1)*exp(beta(2)./x); --一定要./
或function f=zhang1(beta,x)
a=beta(1);
b=beta(2);
f=a*exp(b./x);-------------一定要点除
2、输入数据:
>> x=2:16;
>> y=[6.42 8.20 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76];
>> beta0=[8 2]';----初值[1,1]也可以
3、求回归系数:
>> [beta,r ,J]=nlinfit(x',y','volum',beta0); %beta0初值为列/行向量都可以,还是为列吧。
>> beta
beta =
11.6037
-1.01
即得回归模型为:
4、预测及作图:
>> [YY,delta]=nlpredci('volum',x',beta,r ,J)
>>plot(x,y,'k+',x,YY,'r')
或>> plot(x,y,'ro')
>> hold on
>> xx=2:0.05:16;
>> yy=beta(1)*exp(beta(2)./xx);
>> plot(xx,yy,'g')
又或>> plot(x,y,'ro') >> hold on
>> xx=2:0.05:16;
>> yy=volum(beta,xx);--------通过调用用户自编的函数
>> plot(xx,yy,'g')
>> [beta,r ,J]=nlinfit(x',y','volum',[1,1]); %下面换了多个初值,结果都是一样的。
>> beta
beta =
11.6037 -1.01
>> [beta,r ,J]=nlinfit(x',y','volum',[1,5]);
>> beta
beta =
11.6037 -1.01
>> [beta,r ,J]=nlinfit(x',y','volum',[10,5]);
>> beta
beta =
11.6037 -1.01
>> [beta,r ,J]=nlinfit(x',y','volum',[10,50]);
>> beta
beta =
11.6037 -1.01
以下用来lsqcurvefit求解,结果是一样的。
>> [beta,a,b,exitflag]=lsqcurvefit('volum',[8,2],x',y')
Optimization terminated: relative function value
changing by less than OPTIONS.TolFun.
beta =
11.6037 -1.01
a = 0.8627
b = 0.3959
-0.0615
-0.6868
-0.1208
0.0179
-0.0328
0.2284
0.3197
-0.0577
-0.0563
0.0190
-0.1083
0.1544
-0.0910
0.0970
exitflag =
3
>> [beta,a,b,exitflag]=lsqcurvefit('volum',[1,1],x',y') %换不同的初值,结果是一样的。
beta = 11.6037 -1.01
exitflag = 3
>> [beta,a,b,exitflag]=lsqcurvefit('volum',[10,1],x',y')
beta = 11.6037 -1.01
exitflag = 3
>> [beta,a,b,exitflag]=lsqcurvefit('volum',[10,5],x',y')
beta = 11.6037 -1.01
exitflag = 3
>> [beta,a,b,exitflag]=lsqcurvefit('volum',[10,50],x',y')
beta = 11.6037 -1.01
exitflag = 3
例4 财政收入预测问题:财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据,试构造预测模型。
财政收入预测问题:财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据,试构造预测模型。
年份 | 国民收入(亿元) | 工业总产值(亿元) | 农业总产值(亿元) | 总人口(万人) | 就业人口(万人) | 固定资产投资(亿元) | 财政收入(亿元) |
1952 | 598 | 349 | 461 | 57482 | 20729 | 44 | 184 |
1953 | 586 | 455 | 475 | 58796 | 213 | 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 | 653 | 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 | 9 | 461 | 67295 | 25110 | 66 | 230 |
1963 | 779 | 1046 | 514 | 69172 | 260 | 85 | 266 |
19 | 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 | 17 | 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 | 5 |
1971 | 1780 | 3156 | 790 | 85229 | 35620 | 355 | 638 |
1972 | 1833 | 3365 | 7 | 87177 | 35854 | 354 | 658 |
1973 | 1978 | 3684 | 855 | 211 | 36652 | 374 | 691 |
1974 | 1993 | 3696 | 1 | 90859 | 37369 | 393 | 655 |
1975 | 2121 | 4254 | 932 | 92421 | 38168 | 462 | 692 |
1976 | 2052 | 4309 | 955 | 93717 | 38834 | 443 | 657 |
1977 | 21 | 4925 | 971 | 94974 | 39377 | 454 | 723 |
1978 | 2475 | 5590 | 1058 | 96259 | 39856 | 550 | 922 |
1979 | 2702 | 6065 | 1150 | 97542 | 40581 | 5 | 0 |
1980 | 2791 | 6592 | 1194 | 98705 | 416 | 568 | 826 |
1981 | 2927 | 6862 | 1273 | 100072 | 73280 | 496 | 810 |
y= ax1+bx2+cx3+dx4+ex5+fx6
使用非线性回归方法求解。
1. 对回归模型建立M文件model.m如下:
function yy=model(beta0,X) %一定是两个参数,第一个为系数数组,b(1),b(2),…b(n) %分别代表每个系数,而第二个参数代表所有的自变量, %是一个矩阵,它的每一列分别代表一个自变量。
a=beta0(1);
b=beta0(2); %每个元素
c=beta0(3);
d=beta0(4);
e=beta0(5);
f=beta0(6);
x1=X(:,1); %每一列
x2=X(:,2);
x3=X(:,3);
x4=X(:,4);
x5=X(:,5);
x6=X(:,6);
yy=a*x1+b*x2+c*x3+d*x4+e*x5+f*x6;
2. 主程序liti6.m如下:
X=[598.00, 349.00 ,461.00, 57482.00, 20729.00, 44.00;----可以直接整个从WORD拷贝过来
586, 455, 475, 58796, 213, ;
707, 520, 491, 60266, 21832, 97;
737, 558, 529, 61465, 22328, 98;
825, 715, 556, 62828, 23018, 150;
837, 798, 575, 653, 23711, 139;
1028, 1235, 598, 65994, 26600, 256;
1114, 1681, 509, 67207, 26173, 338;
1079, 1870, 444, 66207, 25880, 380;
757, 1156, 434, 65859, 25590, 138;
677, 9, 461, 67295, 25110, 66;
779, 1046, 514, 69172, 260, 85;
943, 1250, 584, 70499, 27736, 129;
1152, 1581, 632, 72538, 28670, 175;
1322, 1911, 687, 74542, 29805, 212;
1249, 17, 697, 76368, 30814, 156;
1187, 1565, 680, 78534, 31915, 127;
1372, 2101, 688, 80671, 33225, 207;
1638, 2747, 767, 82992, 34432, 312;
1780, 3156, 790, 85229, 35620, 355;
1833, 3365, 7, 87177, 35854, 354;
1978, 3684, 855, 211, 36652, 374;
1993, 3696, 1, 90859, 37369, 393;
2121, 4254, 932, 92421, 38168, 462;
2052, 4309, 955, 93717, 38834, 443;
21, 4925, 971, 94974, 39377, 454;
2475, 5590, 1058, 96259, 39856, 550;
2702, 6065, 1150, 97542, 40581, 5;
2791, 6592, 1194, 98705, 416, 568;
2927, 6862, 1273, 100072, 73280, 496];
y=[184.00 216.00 248.00 254.00 268.00 286.00 357.00 444.00 506.00 ...
271.00 230.00 266.00 323.00 393.00 466.00 352.00 303.00 447.00 ...
5.00 638.00 658.00 691.00 655.00 692.00 657.00 723.00 922.00 ...
0.00 826.00 810.0]';
beta0=[0.50 -0.03 -0.60 0.01 -0.02 0.35];----也可以用5个1
betafit = nlinfit(X,y,'model',beta0)
结果为:
betafit =
0.5243
-0.0294
-0.6304
0.0112
-0.0230
0.3658
(结果也可能是:0.3459 -0.0180 -0.3700 0.0030 -0.0020 0.4728)
即y= 0.5243x1-0.0294x2-0.6304x3+0.0112x4-0.0230x5+0.3658x6
此题也可以用regress来求解(我自己做的,不一定对)----结果有些不同,含有一个常数
>> clear
>> x=xlsread('cz.xls'); %已经把所有的有效数据拷入到cd.xls文件中去了。
>> y=x(:,7);
>> x(:,7)=[ ];
>> z=ones(30,1);
>> x=[z,x];
>> [b,bint,r,rint,states]=regress(y,x);
>> b,states
b =
159.1440
0.4585
-0.0112
-0.5125
0.0008
-0.0028
0.3165
stats =
1.0e+003 *
0.0010 0.2283 0 1.0488
X,y的原始数据:
x = 598 349 461 57482 20729 44
586 455 475 58796 213
707 520 491 60266 21832 97
737 558 529 61465 22328 98
825 715 556 62828 23018 150
837 798 575 653 23711 139
1028 1235 598 65994 26600 256
1114 1681 509 67207 26173 338
1079 1870 444 66207 25880 380
757 1156 434 65859 25590 138
677 9 461 67295 25110 66
779 1046 514 69172 260 85
943 1250 584 70499 27736 129
1152 1581 632 72538 28670 175
1322 1911 687 74542 29805 212
1249 17 697 76368 30814 156
1187 1565 680 78534 31915 127
1372 2101 688 80671 33225 207
1638 2747 767 82992 34432 312
1780 3156 790 85229 35620 355
1833 3365 7 87177 35854 354
1978 3684 855 211 36652 374
1993 3696 1 90859 37369 393
2121 4254 932 92421 38168 462
2052 4309 955 93717 38834 443
21 4925 971 94974 39377 454
2475 5590 1058 96259 39856 550
2702 6065 1150 97542 40581 5
2791 6592 1194 98705 416 568
2927 6862 1273 100072 73280 496
>> y
y =
184
216
248
254
268
286
357
444
506
271
230
266
323
393
466
352
303
447
5
638
658
691
655
692
657
723
922
0
826
810
help nlinfit/ help nlinfit
Examples:
Use @ to specify MODELFUN:
load reaction;
beta = nlinfit(reactants,rate,@mymodel,beta);
where MYMODEL is a MATLAB function such as:
function yhat = mymodel(beta, x)
yhat = (beta(1)*x(:,2) - x(:,3)/beta(5)) ./ ...
(1+beta(2)*x(:,1)+beta(3)*x(:,2)+beta(4)*x(:,3));
Examples
--------
FUN can be specified using @:
nlintool(x, y, @myfun, b0)
where MYFUN is a MATLAB function such as:
function yhat = myfun(beta, x)
b1 = beta(1);
b2 = beta(2);
yhat = 1 ./ (1 + exp(b1 + b2*x));
FUN can also be an inline object:
fun = inline('1 ./ (1 + exp(b(1) + b(2)*x))', 'b', 'x')
nlintool(x, y, fun, b0)
每一列为一个变量,所以有n个自变量,就有n列
X1,x2,…xn, y都为列向量
此即非线性回归或曲线回归问题
(需要配曲线)
配曲线的一般方法是:
(一)先对两个变量x和y 作n次试验观察得画出散点图,
散点图
(二)根据散点图确定须配曲线的类型.
通常选择的六类曲线如下:
(1)双曲线
(2)幂函数曲线y=a, 其中x>0,a>0
(3)指数曲线y=a其中参数a>0.
(4)倒指数曲线y=a其中a>0,
(5)对数曲线y=a+blogx,x>0
(6)S型曲线
(三)然后由n对试验数据确定每一类曲线的未知参数a和b. (采用的方法是通过变量代换把非线性回归化成线性回归,即采用
非线性回归线性化的方法.)
解例2.由散点图我们选配倒指数曲线y=a
根据线性化方法,算得
由此
最后得
作 业
1、考察温度x对产量y的影响,测得下列10组数据:
温度(℃) | 20 | 25 | 30 | 35 | 40 | 45 | 50 | 55 | 60 | 65 |
产量(kg) | 13.2 | 15.1 | 16.4 | 17.1 | 17.9 | 18.7 | 19.6 | 21.2 | 22.5 | 24.3 |
2、某零件上有一段曲线,为了在程序控制机床上加工这一零件,需要求这段曲线的解析表达式,在曲线横坐标xi处测得纵坐标yi共11对数据如下:
xi | 0 | 2 | 4 | 6 | 8 | 10 | 12 | 14 | 16 | 18 | 20 |
yi | 0.6 | 2.0 | 4.4 | 7.5 | 11.8 | 17.1 | 23.3 | 31.2 | 39.6 | 49.7 | 61.7 |
3: 在研究化学动力学反应过程中,建立了一个反应速度和反应物
含量的数学模型,形式为
其中是未知参数,是三种反应物(氢,n戊烷,
异构戊烷)的含量,y是反应速度.今测得一组数据如表4,试由
此确定参数,并给出置信区间.的参考值为
(1,0.05, 0.02, 0.1, 2)
序号 | 反应速度y | 氢x1 | n戊烷x2 | 异构戊烷x3 |
1 | 8.55 | 470 | 300 | 10 |
2 | 3.79 | 285 | 80 | 10 |
3 | 4.82 | 470 | 300 | 120 |
4 | 0.02 | 470 | 80 | 120 |
5 | 2.75 | 470 | 80 | 10 |
6 | 14.39 | 100 | 190 | 10 |
7 | 2.54 | 100 | 80 | 65 |
8 | 4.35 | 470 | 190 | 65 |
9 | 13.00 | 100 | 300 | 54 |
10 | 8.50 | 100 | 300 | 120 |
11 | 0.05 | 100 | 80 | 120 |
12 | 11.32 | 285 | 300 | 10 |
13 | 3.13 | 285 | 190 | 120 |
养护时间x | 2 | 3 | 4 | 5 | 7 | 9 | 12 | 14 | 17 | 21 | 28 | 56 |
抗压强度y | 35 | 42 | 47 | 53 | 59 | 65 | 68 | 73 | 76 | 82 | 86 | 99 |