(完整版)Matlab线性回归(拟合)

Matlab 线性回归（拟合）对于多元线性回归模型：e x x y p p ++++=βββΛ110设变量12,,,p x x x y L 的n 组观测值为12(,,,)1,2,,i i ip i x x x y i n =L L ．记 ⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=np n n p p x x x x x x x x x x ΛΛΛΛΛΛΛΛ212222*********，⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=n y y y y M 21，则⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=p ββββM 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 21M 则目标函数可表达为∑=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 矩阵。