如何用matlab进行非线性曲线拟合

在自然科学，工程技术和经济活动等各领域中，经常需要根据实验观测数据（xi，yi），i=1,2，.....,n研究因变量y与自变量x之间的关系。一般来说，变量之间的关系分为两种，一种是确定性的函数关系，另一种是不确定性关系，也称为相关关系。

回归分析是研究变量之间的相关关系的数学工具，主要解决以下几个方面的问题。

（1）根据变量观测数据确定某些变量之间的定量关系式，即建立回归方程并估计其中的未知参数。估计参数的常用方法是最小二乘法。

（2）对求得的回归方程的可信度进行统计检验。

（3）判断自变量对因变量有无影响。特别的，在许多自变共同影响一个因变量的关系中，需要判断哪些自变量的影响是显著的，哪些自变量的影响是不显著的，将影响显著的自变量选入模型中，剔除影响不显著的变量，通常用逐步回归、向前回归和向后回归等方法。

（4）利用所求的回归方程对某一生成过程进行预测和控制。

1.MATLAB中，有三种回归模型类：LinearModel class（线性回归模型类）、NonLinear class（非线性回归模型类）和GeneralizedLinearModle class（广义线性回归），通过调用类的构造函数可以创建类对象，然后调用类对象的各种方法（例如fit和predict方法）作回归分析。

1.1线性回归模型类

1.1.1线性回归模型类的类方法

对于一元或多元线性回归，MATLAB提供了LinearModel类，用户可以根据自己的观测数据，调用LinearModel类的方法创建一个LinearModel类对象，用来求解型回归模型。LinearModel类的方法如下表所示：

>> mdl=LinearModel %创建一个LinearModel类对象

mdl =

Linear regression model:

y ~ 0

Coefficients:

>> methods(mdl) %查看LinearModel类中的方法

类 LinearModel 的方法:

addTerms disp plotAdded plotInteraction random

anova dwtest plotAdjustedResponse plotResiduals removeTerms

coefCI feval plotDiagnostics plotSlice step

coefTest plot plotEffects predict

1.1.2 线性回归模型类的类属性

LinearModel类对象的属性中包含了模型求解的所有结果，

>> properties(mdl) %查询LinearModel类对象mdl�

类 LinearModel 的属性:

MSE %均方误差（残差）

Robust %稳健回归参数

Residuals %残差

Fitted %拟合值

Diagnostics %回归诊断统计量

RMSE %均方根误差（残差）

Steps %逐步回归相关信息

Formula %回归模型公式

LogLikelihood %对数似然函数值

DFE %误差（残差）的自由度

SSE %误差（残差）平方和

SST %y的总离差平方和

SSR %回归平方和

CoefficientCovariance %系数统计值的协方差矩阵

CoefficientNames %

NumCoefficients

NumEstimatedCoefficients

Coefficients

Rsquared

ModelCriterion

VariableInfo

ObservationInfo

Variables

NumVariables

VariableNames

NumPredictors

PredictorNames

ResponseName

NumObservations

ObservationNames

1.2 非线性回归模型类

1.2.1 非线性回归模型类NonLinearModel类的类方法

MATLAB中提供了NonLinearModel类，用来求解一元或多元非线性回归模型类。NonLinearModel类的类方法如下表：

1.2.2非线性回归模型类NonLinearModel的类属性

NonLinearModel类对象的属性中包含了模型求解的所有结果，用户可以通过查询NonLinearModel类对象的指定属性的属性值得到想要的结果。

>> nlm=NonLinearModel

nlm =

Nonlinear regression model:

y ~ 0

Coefficients:

>> methods(nlm) %查看类方法

类 NonLinearModel 的方法:

coefCI disp plotDiagnostics plotSlice random

coefTest feval plotResiduals predict

>> properties(nlm) %查看NonLinearModel类属性

类 NonLinearModel 的属性:

MSE

Iterative

Robust

Residuals

Fitted

RMSE

Diagnostics

WeightedResiduals

Formula

LogLikelihood

DFE

SSE

SST

SSR

CoefficientCovariance

CoefficientNames

NumCoefficients

NumEstimatedCoefficients

Coefficients

Rsquared

ModelCriterion

VariableInfo

ObservationInfo

Variables

NumVariables

VariableNames

NumPredictors

PredictorNames

ResponseName

NumObservations

ObservationNames

2. 一元线性回归

现有全国31个城市的气候情况观测数据，在“气温.xls”文件中，现在根据这些观测数据研究年平均气温和全年光照时数之间的关系。

2.1 数据的散点图

令x表示年平均气温，y表示全面日照时数。由于x和y均是一维变量，可以先从x和y的散点图上直观地观察他们之间的关系，然后再进一步的分析。

>> ClimateData = xlsread('气温.xls'); %读取数据文件

x = ClimateData(:, 1); %读取文件中的第2列，即年平均气温数据

y = ClimateData(:, 5); %读取第5列，即全年日照时数数据

figure;

plot(x, y, 'k.', 'Markersize', 15); %绘制x和y的散点图

xlabel('年平均气温(x)');

ylabel('全年日照时数(y)');

R = corrcoef(x, y) %计算x和y的线性相关系数矩阵R

R =

1.0000 -0.7095

-0.7095 1.0000

从散点图可以看出，有4组数据有些异常，除了这4组数据外，散点图编写x和y的线性趋势比较明显，可以用直线y=b0+b1*x进行拟合。

从相关系数可以看到，x和y的线性相关系数为-0.7095，表明x和y负相关。

2.2 模型的建立与求解

（1）模型的建立

建立y关于x的一元线性回归模型如下：

yi=b0+b1*x+a

a~N(0,c^2), i=1,2,...n

模型中有四个基本假定：线性假设、误差正态性假设、误差方差齐次性假定、误差独立性假定。

（2）调用LinearModel类的fit方法求解模型

mdl1 = LinearModel.fit(x,y) %模型求解

Linear regression model:

y ~ 1 + x1

Estimated Coefficients:

Estimate SE tStat pValue

________ ______ _______ __________

(Intercept) 3115.4 223.06 13.967 2.0861e-14

x1 -76.962 14.197 -5.4211 7.8739e-06

Number of observations: 31, Error degrees of freedom: 29

Root Mean Squared Error: 383

R-squared: 0.503, Adjusted R-Squared 0.486

F-statistic vs. constant model: 29.4, p-value = 7.87e-06

从输出的结果看，常数项b0和回归系数b1,的估计值分别为3115.4和-76.962，从而可以写出线性回归方程为

y=3115.4-76.962*x

对回归直线进行显著性检验的p=7.8739e-06<0.05,认为是y与x的线性关系是显著的。

（3）调用LinearModel类的plot方法绘制拟合效果图

下面调用LinearModel类的plot方法绘制拟合效果图。

（4）预测

给定自变量x的值，可以调用LinearModel类对象的predict方法计算因变量y的预测值。例如，给定年平均气温x=5,25，计算全年日照时数y的预测值？

xnew = [5,25]'; %定义新的自变量，必须是列向量矩阵

ynew = mdl1.predict(xnew) %计算预测值

ynew =

1.0e+03 *

2.7306

1.1913

2.3 回归诊断

回归诊断主要包括以下内容：

异常点和强影响点诊断，查找数据集中的异常点（离群点）和强影响点，对模型做出改进；

残差分析，用来验证模型的基本假定，包括模型线性诊断，误差正态性诊断，误差方差齐次性诊断和误差独立性诊断

多重共线性诊断，对应多元线性回归，检验自变量之间是否存在共线性关系。

（1）查找异常点和强影响点

数据中的异常点是指远离数据中心的观测点，又称离群点，强影响点是指数据集中对回归方差参数估计结果有较大影响的观测点，通过剔除这些异常点和某些强影响点，可对模型做出改进。 LinearModel类对象的Residuals属性值中列出了标准化残差和学生化残差值，（用来查找异常点）

在回归分析中，测定值与按回归方程预测的值之差称为残差，以δ表示。残差δ遵从正态分布N(0，σ2)。（δ-残差的均值）/残差的标准差，称为标准化残差，以δ*表示。δ*遵从标准正态分布N(0，1)。实验点的标准化残差落在(-2，2)区间以外的概率≤0.05。若某一实验点的标准化残差落在(-2，2)区间以外，可在95%置信度将其判为异常实验点，不参与回归直线拟合

LinearModel类对象的Diagnostics属性值中包含有杠杆值，Cook距离，Covration统计量，Dffits统计量，Dfbeta统计量等回归诊断结果，这些是用来查找强影响点的。可以调用LinearModel类对象的plotDIagnostics方法绘制个统计量对应的回归诊断图，借助回归诊断图可以直观的查找异常点和强影响点。

Res = mdl1.Residuals; %查询残差值

Res_Stu = Res.Studentized; %学生化残差

Res_Stan = Res.Standardized; %标准残差

figure;

subplot(2,3,1);

plot(Res_Stu,'kx'); %绘制学生化残差图

refline(0,-2); %在y=2和y=-2处画直线，超过这个范围即为异常点

refline(0,2);

title('(a) 学生化残差图')

xlabel('观测序号');ylabel('学生化残差');

subplot(2,3,2);

mdl1.plotDiagnostics('cookd'); %绘制Cook距离图，直线以上的点为强影响点

title('(b) Cook距离图')

xlabel('观测序号');ylabel('Cook距离');

subplot(2,3,3);

mdl1.plotDiagnostics('covratio'); %线外的为强影响点

title('(c) Covratio统计量图');

xlabel('观测序号');ylabel('Covratio统计量');

subplot(2,3,4);

plot(Res_Stan,'kx'); %绘制标准化残差图

refline(0,-2); %在y=2和y=-2处画直线，超过这个范围即为异常点

refline(0,2);

title('(d) 标准化残差图');

xlabel('观测序号');ylabel('标准化残差');

subplot(2,3,5);

mdl1.plotDiagnostics('dffits'); %绘制Dffits统计量图，超出点为强影响点

title('(e) Dffits统计量图');

xlabel('观测序号');ylabel('Dffits统计量');

subplot(2,3,6);

mdl1.plotDiagnostics('leverage'); %绘制杠杆值图，超出点为强影响点

title('(f) 杠杆值图');

xlabel('观测序号');ylabel('杠杆值');

（2）模型改进

将检测到的4组异常数据剔除后重新做一元线性回归分析，对模型做出改进。

id = find(abs(Res_Stu)>2); %查找异常值的序号

mdl2 = LinearModel.fit(x,y,'Exclude',id) %去除异常值，重新求解

figure;

mdl2.plot; %绘制拟合效果图

xlabel('年平均气温(x)');

ylabel('全年日照时数(y)');

title('');

legend('剔除异常数据后散点','回归直线','置信区间'); %图例

Linear regression model:

y ~ 1 + x1

Estimated Coefficients:

Estimate SE tStat pValue

________ ______ _______ __________

(Intercept) 2983.8 121.29 24.601 4.8701e-19

x1 -63.628 7.7043 -8.2587 1.3088e-08

Number of observations: 27, Error degrees of freedom: 25

Root Mean Squared Error: 201

R-squared: 0.732, Adjusted R-Squared 0.721

F-statistic vs. constant model: 68.2, p-value = 1.31e-08

剔除异常值后的回归直线方程为y=2983.8-63.628*x

对回归直线进行显著性检验的p值为1.3088e-08,可知y与x的线性关系更为显著，拟合效果图如下：

为了便于比较，做出原始数据散点、原始数据对应的回归线和剔除异常数据后的回归直线图，

figure; %新建一个图形窗口

plot(x,y,'ko'); %画出原始数据散点

hold on; %图形叠加

xnew=sort(x); %将x从小到大排序，为了画图的需要

yhat1=mdl1.predict(xnew); %计算模型1的拟合值

yhat2=mdl2.predict(xnew); %计算模型2的拟合值

plot(xnew,yhat1,'r--','linewidth',3); %画原始数据对应的回归线，红色虚线

plot(xnew,yhat2,'b','LineWidth',3); %画出剔除异常数据后的回归直线，蓝色实线

legend('原始数据散点','原始数据回归直线','剔除异常数据后回归直线'); %图例

xlabel('年平均气温');

ylabel('全年日照时数')

由于受异常数据的影响，两次回归结果并不相同。

（3）残差分析

在回归诊断中，常借助残差图来验证模型的基本假定是否成立。常用的残差图包括残差序列图，残差与拟合值图、残差直方图、残差正态概率图、残差与滞后残差图。可以通过调用LinearModel类对象的plotResiduals方法绘制各种残差图。

figure;

subplot(2,3,1);

mdl2.plotResiduals('caseorder'); %绘制残差值序列图

title('(a) 残差值序列图');

xlabel('观测序号');ylabel('残差');

subplot(2,3,2);

mdl2.plotResiduals('fitted'); %绘制残差与拟合值图

title('(b) 残差与拟合值图');

xlabel('拟合值');ylabel('残差');

subplot(2,3,3);

plot(x,mdl2.Re;kx'); %绘制残差与自变量图

refline(0,0);

title('(c) 残差与自变量图');

xlabel('自变量值');ylabel('残差');

subplot(2,3,4);

mdl2.plotResiduals('histogram'); %绘制残差直方图

title('(d) 残差直方图');

xlabel('残差r');ylabel('f(r)');

subplot(2,3,5);

mdl2.plotResiduals('probability'); %绘制残差正态概率图

title('(e) 残差正态概率图');

xlabel('残差');ylabel('概率');

subplot(2,3,6);

mdl2.plotResiduals('lagged'); %绘制操作与滞后残差图

title('(f) 残差与滞后残差图');

xlabel('滞后残差');ylabel('残差');

（a）残差值序列图，横坐标为观测序号，纵坐标为残差值，可以看出个观测对应的残差值随机地在水平轴y=0上下无规则地波动，说明残差值间是相互独立的。如果残差的分布有一定的规律，则说明残差间不独立。

（b）残差与拟合值图，横坐标为拟合值，纵坐标为残差值，可以看出残差基本分布在左右等宽的水平条带内，说明残差值是等方差。如果残差分布呈现喇叭口形，则说明残差不满足方差齐性假定，此时应对因变量y做某种变换（如取平方根、取对数、去倒数等），然后重新拟合。

（c）残差与自变量图，横坐标为自变量的值，纵坐标为残差值，可以看出残差基本分布在左右等宽的水平条带内，说明线性模拟与数据拟合较好。

（d）残差直方图，残差直方图反映了残差的分布。因为这里数据太少，不能根据残差直方图验证残差的正态性。

（e）残差正态概率图，用来检验残差是否服从正态分布。从图中可以看出残差基本服从正态分布

（f）残差与滞后残差图，横坐标为滞后残差，纵坐标为残差，用来检验残差间是否存在自相关性。从图中可以看出散点均匀分布在四个象限内，说明残差间不存在自相关性。

2.4 稳健回归

默认情况下，fit函数的‘RobustOpts’参数值为‘off’，此时fit函数利用普通最小二乘法估计模型中的参数，参数的估计值受异常值的影响比较大。若将fit函数的‘RobustOpt’参数值设置为‘on’，则可采用加权最小二乘法估计模型中的参数，结果受异常值的影响就比较小。

下面是稳健回归的MATLAB实现

mdl3 = LinearModel.fit(x,y,'RobustOpts','on')

mdl3 =

Linear regression model (robust fit):

y ~ 1 + x1

Estimated Coefficients:

Estimate SE tStat pValue

________ ______ ______ __________

(Intercept) 3034.8 182.01 16.674 2.1276e-16

x1 -68.3 11.584 -5.896 2.1194e-06

Number of observations: 31, Error degrees of freedom: 29

Root Mean Squared Error: 313

R-squared: 0.551, Adjusted R-Squared 0.535

F-statistic vs. constant model: 35.5, p-value = 1.78e-06

稳健回归得出的回归方差为y=3034.8-68.3*x。

下面对比非稳健拟合和稳健拟合，可以看出加权最小二乘拟合的稳健性。

mdl3 = LinearModel.fit(x,y,'RobustOpts','on')

xnew = sort(x); %将x从小到大排序，为了后面画图的需要

yhat1 = mdl1.predict(xnew); %计算拟合值（非稳健拟合）

yhat3 = mdl3.predict(xnew); %计算拟合值（稳健拟合）

plot(x, y, 'ko');

hold on;

plot(xnew, yhat1, 'r--','linewidth',3);

plot(xnew, yhat3, 'linewidth', 3);

legend('原始数据散点','非稳健拟合回归直线','稳健拟合回归直线');

xlabel('年平均气温(x)');

ylabel('全年日照时数(y)');

非稳健拟合是基于普通最小二乘拟合，而稳健拟合是基于加权最小二乘拟合。从图中可以看出，通过加权可以消除异常值的影响，增强拟合的文件性。