非线性模型参数估计

1常见模型

1.1指数模型估计

指数模型数学表达式：

两边取对数，化为线性形式：

%%非线性回归，指数模型
x=[0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5];
y=[26300,25100,19900,15500,11500,9800,5200,4600,3200,2300,1700,1200,900,700,600,500];
plot(x,y,'r.');  %%绘制散点图
xlabel('Distance');
ylabel('Average density');
hold on
X=[ones(length(y),1),x'];  %%自变量矩阵
Y=log(y'); %%因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);  %%回归分析
R2=Stats(1);  %%拟合优度
a=exp(B(1));  %%模型
b=-B(2);  %%回归系数
f=a*exp(-b*x);  %%模型表达式
plot(x,f,'b-');  %%添加趋势线
hold off
s=sqrt(sumsqr(y-f)/(length(f)-2));  %%计算标准误差
a,b,R2,s  %%输出主要变量的计算结果

 

指数模型拟合

%%非线性回归，指数模型
x=[0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5];
y=[26300,25100,19900,15500,11500,9800,5200,4600,3200,2300,1700,1200,900,700,600,500];
plot(x,y,'r.');  %%绘制散点图
xlabel('Distance');
ylabel('Average density');
hold on
beta0=[0,0];  %%设定迭代的初始值
Option=statset('MaxIter',200);  %%设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,Option); %%非线性拟合
a=B(1);  %%模型常数
b=-B(2); %%回归系数
f=a*exp(-b*x);  %%模型表达式
plot(x,f,'b-');
hold off
s=sqrt(sumsqr(y-f)/(length(f)-2));  %%计算标准误差
a,b,s  %%输出主要变量结果

创建myfun.m文件如下，然后将其保存到工作路径下

%%指数模型的拟合
%%myfun.m
function yhat = myfun(beta,x)
b1=beta(1);
b2=beta(2);
yhat = b1*exp(b2*x);
%%以上模型通过编辑窗口保存在Matlab的worck文件夹中

执行结果如下：

  

1.2对数模型

对数模型的数学表达式：

%%对数模型的线性回归
x=[100.6,103.5,231.3,120.4,230.4,233.3,162.7,236.3,158.7,145.2,207.0,203.3,433.5,372.9,525.3,629.2,963.4,608.8,876.7,832.0,...
    703.1,872.6,2196.2,2422.4,2230.5,1117.1,2558.6,1190.0,1750.2,3710.0,6050.0,5760.0,4460.0,6618.2,6272.8,3840.0,6926.3,3580.0,5817.2,6610.0];
y=[2.6,4.0,6.7,8.9,10.2,12.6,13.6,18.0,21.0,22.0,26.1,28.0,31.0,32.0,34.0,36.4,38.6,40.8,43.0,43.5,...
    46.4,48.0,50.0,52.7,53.3,59.0,60.1,62.0,63.4,67.0,69.0,70.0,72.0,74.0,76.0,78.0,80.7,82.0,86.4,88.0];
plot(x,y,'r.');  %%绘制散点图
xlabel('Per capita income');  %%横轴标签
ylabel('Percent Urban');  %%纵轴标签
hold on
X=[ones(length(y),1),log(x')];   %%自变量矩阵
Y=y';   %%因变量矩阵
[B,Bint,E,Eint,Stats]=regress(Y,X);  %%回归分析
R2=Stats(1);   %%拟合优度
a=B(1);   %%模型常数
b=B(2);  %%回归系数
f=a+b*log(x);  %%模型表达式
s=sqrt(sumsqr(y-f)/(length(f)-2));   %%计算标准误差
f=a+b*log(sort(x));    %%模型表达式
plot(sort(x),f,'b-');   %%添加趋势线
hold off
a,b,R2,s   %%输出主要结果
x=sort(x);  %%自变量排序

执行结果：

    

对数模型的非线性拟合

%%对数模型的非线性拟合

x=[100.6,103.5,231.3,120.4,230.4,233.3,162.7,236.3,158.7,145.2,207.0,203.3,433.5,372.9,525.3,629.2,963.4,608.8,876.7,832.0,...
    703.1,872.6,2196.2,2422.4,2230.5,1117.1,2558.6,1190.0,1750.2,3710.0,6050.0,5760.0,4460.0,6618.2,6272.8,3840.0,6926.3,3580.0,5817.2,6610.0];
y=[2.6,4.0,6.7,8.9,10.2,12.6,13.6,18.0,21.0,22.0,26.1,28.0,31.0,32.0,34.0,36.4,38.6,40.8,43.0,43.5,...
    46.4,48.0,50.0,52.7,53.3,59.0,60.1,62.0,63.4,67.0,69.0,70.0,72.0,74.0,76.0,78.0,80.7,82.0,86.4,88.0];
plot(x,y,'r.');  %%绘制散点图
xlabel('Per capita income');  %%横轴标签
ylabel('Percent Urban');  %%纵轴标签
hold on
beta0=[0 0];  %%设定迭代初始值
O=statset('MaxIter',200);   %%设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);  %%非线性拟合
a=B(1);   %%模型常数
b=B(2);  %%回归系数
f=a+b*log(x);  %%模型表达式
s=sqrt(sumsqr(y-f)/(length(f)-2));   %%计算标准误差
f=a+b*log(sort(x));    %%模型表达式
plot(sort(x),f,'b-');   %%添加趋势线
hold off
a,b,s  %%输出主要结果

%%myfun.m
function yhat=myfun(beta,x)
b1=beta(1);
b2=beta(2);
yhat=b1+b2*log(x);
%%以上函数模型另存为.m文件

执行结果：

      

1.3 幂指数模型

指数模型数学表达式：

x=[3.48,3.69,1.43,2.05,3.05,3.19,2.43,2.40,2.72,2.99,3.78,1.33,1.67,3.14,2.04,1.77,0.59,0.69,0.50,0.69,0.63,0.58,0.86,0.41,1.23];
y=[38.83,43.92,9.14,16.66,36.16,38.66,17.74,19.46,23.00,29.75,51.19,6.60,9.04,33.27,17.61,13.37,2.04,2.22,1.46,1.92,1.86,1.69,3.31,1.13,6.74];
plot(x,y,'r.');  %%绘制散点图
xlabel('Perimeter');  %%横轴标签
ylabel('Area');  %%纵轴标签
hold on
X=[ones(length(y),1),log(x')];  %%自变量矩阵
Y=log(y');   %%因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);   %%回归分析
R2=Stats(1);  %%拟合优度
a=exp(B(1));   %%模型常数
b=B(2);  %%回归系数
f=a*x.^b;  %%模型表达式
s=sqrt(sumsqr(y-f)/(length(f)-2));   %%计算标准误差
f=a*(sort(x)).^b;    %%模型表达式
plot(sort(x),f,'b-');   %%添加趋势线
hold off
a,b,R2,s  %%输出主要结果

执行结果：

      

幂指数模型的非线性拟合

x=[3.48,3.69,1.43,2.05,3.05,3.19,2.43,2.40,2.72,2.99,3.78,1.33,1.67,3.14,2.04,1.77,0.59,0.69,0.50,0.69,0.63,0.58,0.86,0.41,1.23];
y=[38.83,43.92,9.14,16.66,36.16,38.66,17.74,19.46,23.00,29.75,51.19,6.60,9.04,33.27,17.61,13.37,2.04,2.22,1.46,1.92,1.86,1.69,3.31,1.13,6.74];
plot(x,y,'r.');  %%绘制散点图
xlabel('Perimeter');  %%横轴标签
ylabel('Area');  %%纵轴标签
hold on
beta0=[0 0];   %%设定迭代初始值
O=statset('MaxIter',200);  %%设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);  %%非线性拟合
a=B(1);   %%模型常数
b=B(2);  %%回归系数
f=a*x.^b;  %%模型表达式
s=sqrt(sumsqr(y-f)/(length(f)-2));   %%计算标准误差
f=a*(sort(x)).^b;    %%模型表达式
plot(sort(x),f,'b-');   %%添加趋势线
hold off
a,b,s  %%输出主要结果

%%myfun.m
function yhat=myfun(beta,x)
b1=beta(1);
b2=beta(2);
yhat=b1*x.^b2;
%%以上函数模型另存为.m文件

执行结果：

     

1.4 多项式模型

n阶多项式函数形式如下：

x=[3.4,1.8,4.6,2.3,3.1,5.5,0.7,3.0,2.6,4.3,2.1,1.1,6.1,4.8,3.8];
y=[26.2,17.8,31.3,23.1,27.5,36.0,14.1,22.3,19.6,31.3,24.0,17.3,43.2,36.4,26.1];
plot(x,y,'r.');  %%绘制散点图
xlabel('Distance');  %%横轴标签
ylabel('Loss');  %%纵轴标签
hold on
X=[ones(length(y),1),x',x'.^2];   %%自变量矩阵
Y=y';   %%因变量矩阵
[B,Bint,E,Eint,Stats]=regress(Y,X);  %%回归分析
R2=Stats(1);    %%拟合优度
f=B(1)+B(2)*x+B(3)*x.^2;  %%模型表达式
s=sqrt(sumsqr(y-f)/(length(f)-3));   %%计算标准误差
x=sort(x);    %%自变量排序
f=B(1)+B(2)*x+B(3)*x.^2;  %%模型表达式
plot(x,f,'b-');   %%添加趋势线
hold off
B,R2,s  %%输出主要结果

执行结果：

      

多项式拟合：

x=[3.4,1.8,4.6,2.3,3.1,5.5,0.7,3.0,2.6,4.3,2.1,1.1,6.1,4.8,3.8];
y=[26.2,17.8,31.3,23.1,27.5,36.0,14.1,22.3,19.6,31.3,24.0,17.3,43.2,36.4,26.1];
plot(x,y,'r.');  %%绘制散点图
xlabel('Distance');  %%横轴标签
ylabel('Loss');  %%纵轴标签
hold on
N=2;   %%确定多项式拟合阶数
[P,S]=polyfit(x,y,N);   %%N阶多项式拟合
lx=linspace(min(x),max(x));  %%限定趋势线的长度范围
z=polyval(P,lx);    %%计算分割点处多项式的值
plot(lx,z,'b-');   %%添加趋势线
hold off
e=y-polyval(P,x);  %%计算残差
s=sqrt(sumsqr(e)/(length(y)-N-1));   %%计算标准误差
P,s,e  %%输出主要结果

执行结果：

 

参考文献：陈彦光.2012.基于Matlab的地理数据分析.高等教育出版社