Matlab 线性回归(拟合)对于多元线性回归模型:e x x y p p ++++=βββ 110设变量12,,,p x x x y 的n 组观测值为12(,,,)1,2,,i i ip i x x x y i n =.记 ⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=np n n p p x x x x x x x x x x 212222111211111,⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=n y y y y 21,则⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=p ββββ 10 的估计值为 y x x x b ')'(ˆ1-==β(11.2) 在Matlab 中,用regress 函数进行多元线性回归分析,应用方法如下:语法:b = regress(y, x)[b, bint, r, rint, stats] = regress(y, x)[b, bint, r, rint, stats] = regress(y, x, alpha)b = regress(y, x),得到的1+p 维列向量b 即为(11.2)式给出的回归系数β的估计值.[b, bint, r, rint, stats]=regress(y, x) 给出回归系数β的估计值b ,β的95%置信区间((1)2p +⨯向量)bint ,残差r 以及每个残差的95%置信区间(2⨯n 向量)rint ;向量stats 给出回归的R 2统计量和F 以及临界概率p 的值.如果i β的置信区间(bint 的第1i +行)不包含0,则在显著水平为α时拒绝0i β=的假设,认为变量i x 是显著的.[b, bint, r, rint, stats]=regress(y, x, alpha) 给出了bint 和rint 的100(1-alpha)%的置信区间.三次样条插值函数的MATLAB 程序matlab 的splinex = 0:10; y = sin(x); %插值点xx = 0:.25:10; %绘图点yy = spline(x,y,xx);plot(x,y,'o',xx,yy)非线性拟合非线性拟合可以用以下命令(同样适用于线形回归分析):1.beta = nlinfit(X,y,fun,beta0)X给定的自变量数据,Y给定的因变量数据,fun要拟合的函数模型(句柄函数或者内联函数形式),beta0函数模型中系数估计初值,beta返回拟合后的系数2.x = lsqcurvefit(fun,x0,xdata,ydata)fun要拟合的目标函数,x0目标函数中的系数估计初值,xdata自变量数据,ydata 函数值数据X拟合返回的系数(拟合结果)nlinfit格式:[beta,r,J]=nlinfit(x,y,’model’, beta0)Beta 估计出的回归系数r 残差J Jacobian矩阵x,y 输入数据x、y分别为n*m矩阵和n维列向量,对一元非线性回归,x为n维列向量。
model是事先用m-文件定义的非线性函数beta0 回归系数的初值例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)x1=x(:,1);x2=x(:,2);x3=x(:,3);yy=beta(1)+beta(2)*x2+beta(3)*x3+beta(4)*(x2.^2)+beta(5)*(x3.^2);(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]';[beta,r,j] = nlinfit(x,y,@myfun,beta0)例题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]);a=beta(1),k1=beta(2),k2=beta(3),m=beta(4)%test the modelxx=min(x):max(x);yy=a+k1*exp(m*xx)+k2*exp(-m*xx);plot(x,y,'o',xx,yy,'r')结果:a = 87.5244k1 = 0.0269k2 = -63.4591m = 0.1083图形:非线性最小二乘(非线性数据拟合)的标准形式为L )x (f )x (f )x (f )x (f min 2m 2221x ++++=其中:L 为常数在MA TLAB5.x 中,用函数leastsq 解决这类问题,在6.0版中使用函数lsqnonlin 。
设⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)x (f )x (f )x (f )x (F m 21 则目标函数可表达为∑=i2i 22x )x (f 21)x (F 21min 其中:x 为向量,F(x)为函数向量。
函数 lsqnonlin格式 x = lsqnonlin(fun,x0) %x0为初始解向量;fun 为)x (f i ,i=1,2,…,m ,fun 返回向量值F ,而不是平方和值,平方和隐含在算法中,fun 的定义与前面相同。
x = lsqnonlin(fun,x0,lb,ub) %lb 、ub 定义x 的下界和上界:b u x b l ≤≤。
x = lsqnonlin(fun,x0,lb,ub,options) %options 为指定优化参数,若x 没有界,则lb=[ ],ub=[ ]。
[x,resnorm] = lsqnonlin(…) % resnorm=sum(fun(x).^2),即解x 处目标函数值。
[x,resnorm,residual] = lsqnonlin(…) % residual=fun(x),即解x 处fun 的值。
[x,resnorm,residual,exitflag] = lsqnonlin(…) %exitflag 为终止迭代条件。
[x,resnorm,residual,exitflag,output] = lsqnonlin(…) %output 输出优化信息。
[x,resnorm,residual,exitflag,output,lambda] = lsqnonlin(…) %lambda 为Lagrage乘子。
[x,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(…) %fun 在解x处的Jacobian 矩阵。
例5-17 求下面非线性最小二乘问题∑=--+101k 2x k x k )e e k 22(21初始解向量为x0=[0.3,0.4]。
解:先建立函数文件,并保存为myfun.m ,由于lsqnonlin 中的fun 为向量形式而不是平方和形式,因此,myfun 函数应由)x (f i 建立:21x k x k k e e k 22)x (f --+= k=1,2,…,10function F = myfun(x)k = 1:10;F = 2 + 2*k-exp(k*x(1))-exp(k*x(2));然后调用优化程序:x0 = [0.3 0.4];[x,resnorm] = lsqnonlin(@myfun,x0)结果为:Optimization terminated successfully:Norm of the current step is less than OPTIONS.TolXx =0.2578 0.2578resnorm = %求目标函数值非线性曲线拟合是已知输入向量xdata 和输出向量ydata ,并且知道输入与输出的函数关系为ydata=F(x, xdata),但不知道系数向量x 。
今进行曲线拟合,求x 使得下式成立:∑-=-i2i i 22x )ydata )xdata ,x (F (21ydata )xdata ,x (F 21min 在MA TLAB5.x 中,使用函数curvefit 解决这类问题。
函数 lsqcurvefit格式 x = lsqcurvefit(fun,x0,xdata,ydata)x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)[x,resnorm] = lsqcurvefit(…)[x,resnorm,residual] = lsqcurvefit(…)[x,resnorm,residual,exitflag] = lsqcurvefit(…)[x,resnorm,residual,exitflag,output] = lsqcurvefit(…)[x,resnorm,residual,exitflag,output,lambda] = lsqcurvefit(…)[x,resnorm,residual,exitflag,output,lambda,jacobian] =lsqcurvefit(…)参数说明:x0为初始解向量;xdata ,ydata 为满足关系ydata=F(x, xdata)的数据;lb 、ub 为解向量的下界和上界b u x b l ≤≤,若没有指定界,则lb=[ ],ub=[ ]; options 为指定的优化参数;fun 为拟合函数,其定义方式为:x = lsqcurvefit(@myfun,x0,xdata,ydata),其中myfun 已定义为 function F = myfun(x,xdata)F = … % 计算x 处拟合函数值fun 的用法与前面相同;resnorm=sum ((fun(x,xdata)-ydata).^2),即在x 处残差的平方和;residual=fun(x,xdata)-ydata ,即在x 处的残差;exitflag 为终止迭代的条件;output 为输出的优化信息;lambda 为解x 处的Lagrange 乘子;jacobian 为解x 处拟合函数fun 的jacobian 矩阵。