【4.0】数学建模中拟合算法详解|内附清晰图片和详细代码实现

一、前言

与插值问题不同，在拟合问题中不需要曲线一定经过给定的点。拟合问题的目标是寻求一个函数（曲线），使得该曲线在某种准则下与所有的数据点最为接近，即曲线拟合的最好（最小损失函数）

插值和拟合的区别：插值算法中，得到的多项式 f ( x ) f(x) f(x)要经过所有的样本点。但是如果样本点过多，得到的多项式会次数很高，从而导致龙格现象。尽管可以采用分段的方法避免这种现象，大使更多的时候更倾向于得到一个确定的曲线，尽管这条曲线不能经过所有的样本点，但只要保证误差足够小即可，这就是拟合的思想（拟合的结果就是得到一个确定的曲线）。

二、最小二乘法

设样本点为 ( x i , y i ) , i = 1 , 2 , … ， n (x_i,y_i),i=1,2,…，n (xi,yi),i=1,2,…，n,假设拟合曲线为 y = k x + b y=kx+b y=kx+b

利用最小二乘法求定义求拟合曲线的方法为：
y ^ i = k x i + b \hat{y}_i=kx_i+b y^i=kxi+b

其中： k ^ i , b ^ i = a r g min ⁡ k , b ( ∑ i = 1 n ( y i − y ^ i ) 2 ) \text{其中：}\hat{k}_i,\hat{b}_i=\underset{k,b}{arg\min}\left( \sum_{i=1}^n{\left( y_i-\hat{y}_i \right) ^2} \right) 其中：k^i,b^i=k,bargmin(i=1∑n(yi−y^i)2)

三、求解最小二乘法

设样本点为 ( x i , y i ) , i = 1 , 2 , … ， n (x_i,y_i),i=1,2,…，n (xi,yi),i=1,2,…，n,设置拟合曲线为 y = k x + b y=kx+b y=kx+b，令拟合值为 y ^ i = k x i + b \hat{y}_i=kx_i+b y^i=kxi+b

那么：
k ^ i , b ^ i = a r g min ⁡ k , b ( ∑ i = 1 n ( y i − y ^ i ) 2 ) = a r g min ⁡ k , b ( ∑ i = 1 n ( y i − k x i − b ) 2 ) \hat{k}_i,\hat{b}_i=\underset{k,b}{arg\min}\left( \sum_{i=1}^n{\left( y_i-\hat{y}_i \right) ^2} \right) =\underset{k,b}{arg\min}\left( \sum_{i=1}^n{\left( y_i-kx_i-b \right) ^2} \right) k^i,b^i=k,bargmin(i=1∑n(yi−y^i)2)=k,bargmin(i=1∑n(yi−kxi−b)2)
令 L = ( ∑ i = 1 n ( y i − k x i − b ) 2 ) L=\left( \sum_{i=1}^n{\left( y_i-kx_i-b \right) ^2} \right) L=(∑i=1n(yi−kxi−b)2)，要找 k k k， b b b使得L最小，在机器学习中 L L L被称之为损失函数，在回归中也常被称之为残差平方和

求解满足要求的 k k k， b b b步骤如下：

四、MATLAB实现最小二乘法拟合

拟合的原始数据为：

x = 4.2,5.9,2.7,3.8,3.8,5.6,6.9,3.5,3.6,2.9,4.2,6.1,5.5,6.6,2.9,3.3,5.9,6,5.6;

y = 8.4,11.7,4.2,6.1,7.9,10.2,13.2,6.6,6,4.6,8.4,12,10.3,13.3,4.6,6.7,10.8,11.5,9.9

五、评价拟合的好坏

几个相关指标

总体平方和SST： S S T = ∑ i = 1 n ( y i − y ˉ ) 2 SST=\sum_{i=1}^n{\left( y_i-\bar{y} \right)}^2 SST=∑i=1n(yi−yˉ)2

误差平方和SSE： S S E = ∑ i = 1 n ( y i − y ^ i ) 2 SSE=\sum_{i=1}^n{\left( y_i-\hat{y}_i \right)}^2 SSE=∑i=1n(yi−y^i)2

回归平方和SSR： S S R = ∑ i = 1 n ( y ^ i − y ˉ ) 2 SSR=\sum_{i=1}^n{\left( \hat{y}_i-\bar{y} \right)}^2 SSR=∑i=1n(y^i−yˉ)2

可以证明得到： S S T = S S E + S S R SST=SSE+SSR SST=SSE+SSR

拟合优度（可决系数）： 0 ⩽ R 2 = S S R S S T = S S T − S S E S S T = 1 − S S T − S S E S S T ⩽ 1 0 \leqslant \,\,R^2=\frac{SSR}{SST}=\frac{SST-SSE}{SST}=1-\frac{SST-SSE}{SST}\,\,\leqslant \,\,1 0⩽R2=SSTSSR=SSTSST−SSE=1−SSTSST−SSE⩽1

说明： R 2 R^2 R2越接近于1，说明误差平方和月接近0，误差越小说明拟合的越好

注意： R 2 R^2 R2只能用于拟合函数是线性函数时（对参数为线性），拟合结果的评价。如果是比较线性函数和其他函数（例如复杂的指数函数）之间拟合的好坏，直接看 S S E SSE SSE指标即可

六、计算拟合优度的代码

y_hat = k*x+b; % y的拟合值
SSR = sum((y_hat-mean(y)).^2)  % 回归平方和
SSE = sum((y_hat-y).^2) % 误差平方和
SST = sum((y-mean(y)).^2) % 总体平方和
SST-SSE-SSR   % 5.6843e-14  =   5.6843*10^-14   matlab浮点数计算的一个误差
R_2 = SSR / SST

七、MATLAB中一个好用的曲线拟合工具箱

八、利用拟合工具预测美国人口

下表给出而来进2个世界的美国人口统计数据（单位：百万人），请使用下面给定的拟合函数预测后30年美国的人口数

年	1790	1800	1810	1820	1830	1840	1850	1860
人口	3.9	5.3	7.2	9.6	12.9	17.1	13.2	31.4
年	1870	1880	1890	1900	1910	1920	1930	1940
人口	38.6	50.2	62.9	76.0	92.0	106.5	123.2	131.7
年	1950	1960	1970	1980	1990	2000
人口	150.7	179.3	204.0	226.5	251.4	281.4

拟合函数如下：

x ( t ) = x m 1 + ( x m 3.9 − 1 ) e − r ( t − 1790 ) x\left( t \right) =\frac{x_m}{1+\left( \frac{x_m}{3.9}-1 \right) e^{-r\left( t-1790 \right)}} x(t)=1+(3.9xm−1)e−r(t−1790)xm

x m x_m xm和 r r r是两个拟合参数， t t t表示年份， x ( t ) x(t) x(t)表示第 t t t年的人口

代码参考：

clear;clc
year = 1790:10:2000;
population = [3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0,106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4];
plot(year,population,'o')
cftool  % 拟合工具箱
% (1) X data 选择 year
% (2) Y data 选择 population
% (3) 拟合方式选择：Custom Equation (自定义方程)
% (4) 修改下方的方框为：x = f(t) = xm/(1+(xm/3.9-1)*exp(-r*(t-1790)))
% (5) 左边的result一栏最上面显示：Fit computation did not converge:即没有找到收敛解，右边的拟合图形也表明拟合结果不理想
% (6) 点击Fit Options，修改非线性最小二乘估计法拟合的初始值(StartPoint), r修改为0.02，xm修改为500
% 有很多同学有疑惑，初始值为什么要这样设置？我们在未来学习微分方程模型和智能算法的课程时再来给大家介绍这里面蕴含的技巧。
% (7) 此时左边的result一览得到了拟合结果：r = 0.02735, xm = 342.4
% (8) 依次点击拟合工具箱的菜单栏最左边的文件—Generate Code(导出代码到时候可以放在你的论文附录)，可以得到一个未命名的脚本文件
% (9) 在这个打开的脚本中按快捷键Ctrl+S，将这个文件保存到当前文件夹。
% (10) 在现在这个文件中调用这个函数得到参数的拟合值和预测的效果
[fitresult, gof] = createFit(year, population)
t = 2001:2030;
xm = 342.4;
r =  0.02735;
predictions = xm./(1+(xm./3.9-1).*exp(-r.*(t-1790)));  % 计算预测值（注意这里要写成点乘和点除,这样可以保证按照对应元素进行计算）
figure(2)
plot(year,population,'o',t,predictions,'.')  % 绘制预测结果图

拟合结果：

更多有关于拟合算法的经典获奖论文，关注公众号，回复，“拟合算法”，即可免费领取！！！