MATLAB数据拟合: 综合实例

[大谦MATLAB，dqmatlab点com]

有理模型拟合

本例用有理模型拟合实测数据。该数据将铜的热膨胀系数描述为开氏温度的函数。首先从文件hahn1.dat中载入热膨胀数据，该文件由工具箱提供。

code.matlab

load hahn1

现在工作空间包括两个新变量，即temp和thermex：

temp是开氏温度组成的矢量。

thermex是铜的热膨胀系数组成的矢量。

将这两个变量导入到曲线拟合工具中，将数据集命名为CuThermEx。现在求拟合该数据集的最佳有理模型。

\[\]

图1-8 曲线拟合工具界面

介绍库模型时讲过，有理模型用多项式的比表示，即

\[y=\frac{{p}_{1}{x}^{n}+{p}_{2}{x}^{n-1}+...+{p}_{n+1}}{{x}^{m}+{q}_{1}{x}^{m-1}+...+{q}_{m}}\]

\[\]

图1-9 拟合界面的设置

其中，n为分子多项式的阶次，m是分母多项式的阶次。注意，有理方程与数据的物理参数无关，相反，它们提供了可以用于内插和外推的简单而且灵活的经验模型。

输入数据，打开曲线拟合工具对话框，如图1-8所示。从数据散点图的形状可以看出，有理模型的最初形式选择分子分母都是二次多项式比较合适。拟合界面的设置如图1-9所示。

数据点、拟合线和残差如图1-10所示。

\[\]

图1-10 拟合结果

二次多项式拟合的结果明显遗漏了最小和最大的预测值。此外，残差点的分布还看得见明显的线形特征，这说明还存在更好的拟合模型。

选择分子与分母都是三次多项式的方程进行拟合，结果如图1-11所示。

\[\]

图1-11 拟合结果

图1-12中的数值结果显示拟合计算不收敛。

\[\]

图1-12 数值结果

“Results”窗口内的消息说明，增加迭代的最大次数可以改进拟合精度。尽管如此，另外使用一个不同的有理方程也许更好，因为当前拟合包含了几个明显的断点。它们是因为分母为零时函数溢出造成的。

下面运用分子是三次多项式，分母是二次多项式的方程进行拟合，效果如图1-13所示。

拟合效果很好，拟合曲线充分体现了整个数据的分布特点，残差随机分布在0附近，可以选择。

\[\]

图1-13 拟合结果图

用自定义方程进行拟合

可以用“Create Custom Equation”对话框定义自己的方程。可以通过下面两种方式打开该对话框：

从“Curve Fitting”对话框中，依次选择菜单Tools->Custom Equation。

在“Fitting”对话框中，从“Type of fit”列表中选择“Custom Equations”选项，然后单击“New Equation”命令按钮。

“Create Custom Equation”界面上有两个面板，一个用于创建线性自定义方程，另一个用于创建非线性自定义方程。这些面板在下面的实例中进行介绍。

1）线性方程：勒让德多项式拟合

本例用几个自定义线性方程拟合数据。数据已经有了，它是从核反应试验中得来的。

有时候用勒让德多项式将变量表示成角度的函数很有用，即

\[y(x)=\sum_{n=0}^{\infty} {{a}_{n}{p}_{n}(x)}\]

其中，Pn(x)为n阶勒让德多项式，x为cos(θα)，an为拟合系数。

可以用Rodrigue公式生成勒让德多项式：

\[{P}_{n}(x)=\frac{1}{{2}^{n}n!}{\left( \frac{d}{dx} \right)}^{n}({x}^{2}-1{)}^{n}\]

表1-2中给出了前4阶勒让德多项式。

表1-2 前4阶勒让德多项式

阶次n	Pn(x)
0	1
1	x
2	`(1/2)(2x2-1)`
3	`(1/2)(5x3-3x)`
4	`(1/2)(35x4-30x2+3)`

本例的第1步是载入文件carbon12alpha.mat中的数据，该数据由工具箱提供。

load carbon12alpha

现在工作空间中包含2个新变量，即angle和counts：

angle为角度矢量，值用弧度表示，以10度为间隔取10度到240度之间的值。

counts为对应于angle中放射性粒子辐射角的粒子数矢量。

在曲线拟合工具箱中导入这两个变量，将数据集命名为C12Alpha。用只有偶数项的4阶勒让德多项式拟合数据。

\[{y}_{1}(x)={a}_{0}+{a}_{2}\left( \frac{1}{2} \right)(3{x}^{2}-1)+{a}_{4}\left( \frac{1}{8} \right)(35{x}^{4}-30{x}^{2}+3)\]

因为勒让德多项式只与自变量和常数有关，所以可以使用“Create Custom Equation”界面上的“Linear Equations”选项卡。对于y1(x)给定的模型，选项卡如图1-14中设置。注意，因为angle用弧度给出，勒让德各项的参数用cos(θα)给定。

\[\]

图1-14 方程设置

拟合曲线和残差如图1-15所示。从图中可以看出，曲线对数据的拟合程度比较高，残差随机分布，看不出任何随机误差。

\[\]

图1-15 拟合曲线和残差

数值拟合结果如图1-16所示。95%置信边界显示与P0(x)和P4(x)有关的系数相当精确，但P2(x)项的系数有较大的不确定性。

\[\]

图1-16 数值拟合结果

为了确定指定偶数项的4阶勒让德多项式能最好地描述数据特征，下面用既有偶数项又有奇数项的多项式进行拟合。

\[{y}_{2}(x)={y}_{1}(x)+{a}_{1}x+{a}_{3}\left( \frac{1}{2} \right)(5{x}^{3}-3x)\]

图1-17所示的数值结果显示，奇数项并没有显著提高拟合精度，而偶数项与前面相比没有多大的变化。这进一步说明了开始的模型是最好的模型。

\[\]

图1-17 数值结果

2) 非线性方程：用傅立叶级数拟合

本例用几个自定义的非线性方程拟合ENSO数据。ENSO数据由复活岛和澳大利亚达尔文之间的月平均大气压差组成。该大气压差造成了南半球信风的形成。ENSO数据具有明显的周期性，可以用傅立叶级数进行描述，即

\[y(x)={a}_{0}+\sum_{i=1}^{\infty} {{a}_{i}\operatorname{cos}{\left( 2\pi\frac{x}{{c}_{i}} \right)}+{b}_{i}\operatorname{sin}{\left( 2\pi\frac{x}{{c}_{i}} \right)}}\]

其中，ai和bi为振幅，ci为周期。本例要解决的问题是数据中存在多少周期。首先假设12个月为一周期，用一个正弦项和一个余弦项拟合数据。

\[{y}_{1}(x)={a}_{0}+{a}_{1}\operatorname{cos}{\left( 2\pi\frac{x}{{c}_{1}} \right)}+{b}_{1}\operatorname{sin}{\left( 2\pi\frac{x}{{c}_{1}} \right)}\]

如果拟合效果不好，添加带唯一周期系数的其他正弦和余弦项，直到拟合效果满意为止。因为三角函数参数中包含了一个未知系数c1，所以方程是非线性的。所以必须用“Create Custom Equation”界面中的“Genteral Equations”选项卡指定方程。图1-18显示了用y1(x)给定的方程。

假设12个月为一周期默认时，系数没有边界限制，初值为0和1之间的随机数指定一个有意义的方程名称假设12个月为一周期默认时，系数没有边界限制，初值为0和1之间的随机数指定一个有意义的方程名称 Inline Image

图1-18 在选项卡中指定方程

图1-19所示的数值结果说明该模型的拟合效果不够好。特别是c1的拟合值小得不合理。因为初值是随机选取的，每次拟合的结果可能会有不同。

\[\]

图1-19 数值结果图1-20 “Fit Options”对话框

为了改进拟合效果，将c1限制在10和14之间。为了定义未知系数的约束，使用“Fit Options”对话框，在“Fitting”对话框中单击“Fit options”命令按钮可以打开它。如图1-20所示。

拟合曲线、残差和数值结果如图1-21所示。

一个周期的拟合数值结果显示12个月为一周期残差说明，至少有1个以上的周期一个周期的拟合数值结果显示12个月为一周期残差说明，至少有1个以上的周期 Inline Image

图1-21 拟合曲线、残差和数值结果

从拟合结果看来，某些点的拟合效果比较好，但很明显，整个数据集的拟合效果不是很好。正如开始所估计的，数值结果说明周期约为12个月。但残差的系统周期分布说明有其他周期应该包含到拟合方程中来。所以，作为第2次尝试，在y1(x)中另外添加一个正弦项和一个余弦项。

\[{y}_{2}(x)={y}_{1}(x)+{a}_{2}\operatorname{cos}{\left( 2\pi\frac{x}{{c}_{2}} \right)}+{b}_{2}\operatorname{sin}{\left( 2\pi\frac{x}{{c}_{2}} \right)}\]

将c2的上下界设置为c1上下界的两倍。拟合结果、残差和数值结果如图1-22所示。

两个周期的拟合数值结果显示，另外一个周期为12个月残差说明，至少有1个以上的周期两个周期的拟合数值结果显示，另外一个周期为12个月残差说明，至少有1个以上的周期 Inline Image

图1-22 添加另一个周期后的拟合结果

拟合结果看起来对大部分数据点都合理，但残差说明应该还包含一个周期到拟合方程中。所以，作为第3次尝试，再添加一个正弦项和余弦项到y2(x)中。

\[{y}_{3}(x)={y}_{2}(x)+{a}_{3}\operatorname{cos}{\left( 2\pi\frac{x}{{c}_{3}} \right)}+{b}_{3}\operatorname{sin}{\left( 2\pi\frac{x}{{c}_{3}} \right)}\]

将c3的下界设置为c1值的3倍。

现在拟合曲线、残差和数值结果如图1-23所示。

三个周期的拟合数值结果显示，三个周期分别为12个月、22个月和44个月对于大部分数据而言，残差分布看起来相当随机三个周期的拟合数值结果显示，三个周期分别为12个月、22个月和44个月对于大部分数据而言，残差分布看起来相当随机 Inline Image

图1-23 添加第3个周期以后的拟合结果

拟合结果说明，现在可以解释ENSO数据中的大部分周期了。尽管残差指示仍然有其他周期存在，但对于大部分数据而言残差还是随机的。也可以改进拟合的振幅值。

总之，数据的傅立叶分析揭示了三个重要的周期。年周期最强，还存在22个月和44个月的周期。

3) 非线性方程：指数背景的高斯模型拟合

本例用非线性自定义模型拟合具有指数衰减特征的数据的两个不太明显的峰值。首先，从文件gauss3.mat载入数据，该数据由工具箱提供。

code.matlab

load gauss3

现在工作空间包含两个新变量，即xpeak和ypeak：

xpeak是由自变量的值组成的矢量。

ypeak是由因变量的值组成的矢量。

将这两个变量导入到曲线拟合工具箱，并接受默认的数据集名称ypeak和xpeak。使用下面的方程拟合数据。

\[y(x)=a{e}^{-bx}+{a}_{1}{e}^{-{\left( \frac{x-{b}_{1}}{{c}_{1}} \right)}^{2}}+{a}_{2}{e}^{-{\left( \frac{x-{b}_{2}}{{c}_{2}} \right)}^{2}}\]

其中，a1为峰的振幅，b1为峰的中心，c1为峰的宽度。因为指数函数参数中包含有未知系数，该方程是非线性的。所以，必须用“Create Custom Equation”对话框的“General Equations”选项卡指定方程，如图1-24所示。

\[\]

图1-24 在“General Equations”选项卡中指定方程

拟合曲线、残差和数值结果如图1-25所示。很明显，拟合结果比较糟糕。

\[\]

图1-25 拟合结果

因为初值是随机选择的，你的拟合结果可能与下面的不同。结果包含下面的警告信息：

Inline Image 图1-26 设置未知系数的初值和约束图1-26 设置未知系数的初值和约束Fit computation did not converge:

Maximum number of function evaluations exceeded. Increasing

MaxFunEvals (in fit options) may allow for a better fit, or

the current equation may not be a good model for the data.

说明计算不收敛，超过了函数计算的最大次数。增大拟合选项MaxFunEvals的值可以获得更好的拟合结果，或者当前方程不是拟合数据的最佳模型。

给系数指定合理的初值点，可以改进本例的拟合效果。对于当前模型，推测初值很容易，因为高斯系数的解释很直接，并且数据的指数特征很明显。另外，因为峰的振幅和宽度不能为负，将a1,a2,c1和c2约束为大于0。

用“Fit Options”对话框定义未知系数的初值和约束，可以通过单击“Fit options”命令按钮打开它。初值和约束如图1-26所示。

数据点、拟合曲线、残差和数值结果如图1-27所示。拟合效果相当不错。

\[\]

图1-27 重新指定初值以后的拟合结果

本页目录

有理模型拟合

用自定义方程进行拟合