Matlab线性回归（拟合）-应用

Matlab线性回归（拟合）

对于多元线性回归模型：

y??0??1x1????pxp?e

设变量

x1,x2,xp,y的n组观测值为

(xi1,xi2,xip,yi)i?1,2,,n．

?1??1记 x?????1?x11x21?xn1?x1p??y1????x22?x2p??y2?，y????， ????????y??xn2?xnp??n?x12??0?????1?则???? 的估计值为

???????p???(x'x)?1x'yb??

在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)，得到的p+1维列向量b即为(11.2)式给出的回归系数β的估计值．

[b, bint, r, rint, stats]=regress(y, x) 给出回归系数β的估计值b，β的95％置信区间（(p+1)*2向量）bint，残差r以及每个残差的95％置信区间（n?2向量）rint；向量stats给出回归的R2统计量和F以及临界概率p的值．

如果

?i的置信区间（bint的第i+1行）不包含0，则在显著水平为?时拒绝

?i?0的假设，认为变量xi是显著的．

[b, bint, r, rint, stats]=regress(y, x, alpha) 给出了bint和rint的100(1-alpha)%的置信区间．

1.三次样条插值函数的MATLAB程序

matlab的spline

x = 0:10; y = sin(x); %插值点 xx = 0:.25:10; %绘图点 yy = spline(x,y,xx);

plot(x,y,'o',xx,yy)

2.非线性拟合

非线性拟合可以用以下命令(同样适用于线性回归分析)： beta = nlinfit(x,y,fun,beta0) x:给定的自变量数据, y:给定的因变量数据,

fun:要拟合的函数模型(句柄函数或者内联函数形式), beta0:函数模型中系数估计初值, beta返回拟合后的系数

x = lsqcurvefit(fun,x0,xdata,ydata) fun要拟合的目标函数,

x0:目标函数中的系数估计初值, xdata:自变量数据, ydata:函数值数据,

x:拟合返回的系数(拟合结果), 2.1 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 model xx=min(x):max(x);

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 图形：

2.2 lsqnonlin

非线性最小二乘（非线性数据拟合）的标准形式为

minxf(x)?f1(x)2?f2(x)2???fm(x)2?L

其中：L为常数

在MATLAB5.x中，用函数leastsq解决这类问题，在6.0版中使用函数lsqnonlin。

?f1(x)??f(x)?2?F(x)???????fm(x)??设

则目标函数可表达为

其中：x为向量，F(x)为函数向量。 函数lsqnonlin

minx1F(x)222?1?fi(x)22i

格式 x = lsqnonlin(fun,x0) %x0为初始解向量；fun为fi(x)，i=1,2,…,m，fun返回向量值F，而不是平方和值，平方和隐含在算法中，fun的定义与前面相同。

x = lsqnonlin(fun,x0,lb,ub) %lb、ub定义x的下界和上界：lb?x?ub。 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求下面非线性最小二乘问题k?1?(2?2k?ekx101?ekx2)2初始解向量为

x0=[0.3, 0.4]。

解：先建立函数文件，并保存为myfun.m，由于lsqnonlin中的fun为向量形式而不是平方和形式，因此，myfun函数应由fi(x)建立：

fk(x)?2?2k?ekx1?ekx2 k=1,2,…,10

function 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.TolX x =

0.2578 0.2578

resnorm = %求目标函数值 lsqcurvefit

非线性曲线拟合是已知输入向量xdata和输出向量ydata，并且知道输入与输出的函数关系为ydata=F(x, xdata)，但不知道系数向量x。今进行曲线拟合，求x使得下式成立：

minx1F(x,xdata)?ydata222?1?(F(x,xdatai)?ydatai)22i

在MATLAB5.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为解向量的下界和上界lb?x?ub，若没有指定界，则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矩阵。 例5-16 求解如下最小二乘非线性拟合问题

已知输入向量xdata和输出向量ydata，且长度都是n，拟合函数为

ydata(i)?x(1)?xdata(i)2?x(2)?sin(xdata(i))?x(3)?xdata(i)3

即目标函数为

minx1(F(x,xdata)?ydata)2ii2?i?1

n23其中：F(x,xdata)?x(1)?xdata?x(2)?sin(xdata)?x(3)?xdata

初始解向量为x0=[0.3, 0.4, 0.1]。

解：先建立拟合函数文件，并保存为myfun.m function F = myfun(x,xdata)

F = x(1)*xdata.^2 + x(2)*sin(xdata) + x(3)*xdata.^3; 然后给出数据xdata和ydata

>>xdata = [3.6 7.7 9.3 4.1 8.6 2.8 1.3 7.9 10.0 5.4];

>>ydata = [16.5 150.6 263.1 24.7 208.5 9.9 2.7 163.9 325.0 54.3]; >>x0 = [10, 10, 10]; %初始估计值

>>[x,resnorm] = lsqcurvefit(@myfun,x0,xdata,ydata) 结果为：

Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun x =

0.2269 0.3385 0.3021 resnorm = 6.2950

进行非线性回归时可使用nlinfit指令，其语法如下：

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

[...] = nlinfit(X, y, fun, beta0, options)

[x,resnorm]= lsqcurvefit(fun,x0,xdata,ydata)； 参数解释：

input：fun——编程者需要拟合的函数 x0——函数系数的初始猜测值 xdata——x坐标的值 ydata——y左边的值

output：x——经拟合的系数

resnorm——the value of the squared 2-norm of the residual sum((fun(x,xdata)-ydata).^2).

example： function fitfunc

xdata=[3.5 7.7 9.3 4.5 8.6 2.8 1.3 7.9 10.1 5.4]; %定义自变量 ydata=[16.5 149.6 263.1 24.6 208.5 9.9 2.7 162.9 322.0 52.3]; %定义因变量

x0=[7,7,7]; %初始估计值

[x,renorm]=lsqcurvefit(@myfun,x0,xdata,ydata); %确定待定系数 disp(x);

disp(renorm);

at x:

本文来源：https://www.bwwdw.com/article/mce2.html