【数学建模笔记 09-2】数学建模的拟合

09-2. 拟合

定义

拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最小。

最小二乘拟合

已知一组点 $(x_i,y_i),i=1,2,\dots ,n$​，$x_i$​ 互不相同，寻求一个函数 $y=f(x)$​，使 $f(x)$ 在某种准则下与所有数据点最为接近。

称
$${\delta }_i=f(x_i)-y_i,i=1,2,\dots ,n$$
为拟合函数 $f(x)$ 在 $x_i$​​ 点的偏差，可以采用偏差平方和最小作为判定准则，即让
$$J=\sum _{i=1}^n(f(x_i)-y_i{)}^2$$
达到最小值。这一原则称最小二乘原则。

线性拟合

令
$$f(x)=a_1{r}_1(x)+a_2{r}_2(x)+\cdots +a_m{r}_m(x),$$
其中 $r_k(x)$ 是事先选定的一组线性无关的函数，$a_k$ 是待定系数。

记
$$J(a_1,\dots ,a_m)=\sum _{i=1}^n[f(x_i)-y_i{]}^2,$$
为求 $a_1,\dots ,a_m$ 使 $J$ 最小，令
$$\frac{\partial J}{\partial a_k}=0,k=1,\dots ,m,$$
得线性方程组
$$\sum _{i=1}^n{r}_j(x_i)[\sum _{k=1}^m{a}_k{r}_k(x_i)-y_i]=0,$$
即
$$\sum _{k=1}^m{a}_k[\sum _{i=1}^n{r}_j(x_i)r_k(x_i)]=\sum _{i=1}^n{r}_j(x_i)y_i,$$

$$j=1,\dots ,m.$$

记
$$R=\left(\begin{array}{ccc}{r}_1(x_1)& \dots & r_m(x_1)\\ ⋮& \ddots & ⋮\\ r_1(x_n)& \dots & r_m(x_n)\end{array}\right),$$

$$A=(a_1,\dots ,a_m{)}^T,Y=(y_1,\dots ,y_n{)}^T.$$

则方程组可表示为
$$R^{T}RA=R^{T}Y,$$
该方程组有唯一解，求解即可。

例子

用 $y=a_0+a_{1}x+a_2{x}^2$ 拟合 $(1,4),(2,10),(3,18),(4,26)$。

有 $r_0(x)=1,r_1(x)=x,r_2(x)=x^2$，故
$$R=\left(\begin{array}{ccc}1& 1& 1\\ 1& 2& 4\\ 1& 3& 9\\ 1& 4& 16\end{array}\right),$$

$$Y=(4,10,18,26{)}^T.$$

由
$$R^{T}RA=R^{T}Y$$
得
$$\left(\begin{array}{ccc}4& 10& 30\\ 10& 30& 100\\ 30& 100& 354\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right)=\left(\begin{array}{c}58\\ 182\\ 622\end{array}\right)$$
求解线性方程组得
$$a_0=-\frac{3}{2},a_1=\frac{49}{10},a_2=\frac{1}{2}.$$
故
$$f(x)=-\frac{3}{2}+\frac{49}{10}x+\frac{1}{2}{x}^2.$$

非线性拟合

对于给定的 $r_k(x)$，如果拟合函数不能以其线性组合的形式出现，则 $f(x)$​ 为非线性函数。对于非线性函数的极小化问题，可用非线性优化方法求解。

Python 代码

使用 numpy 包进行一次、二次、三次多项式拟合上例，代码如下：

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-26
# @ function: 使用 numpy 包进行最小二乘拟合

# %%

import matplotlib.pyplot as plt
import numpy as np

# %%

# 源数据
x = np.array([1, 2, 3, 4])
y = np.array([4, 10, 18, 26])

# 多项式拟合
z1 = np.polyfit(x, y, 1)
z2 = np.polyfit(x, y, 2)
z3 = np.polyfit(x, y, 3)

p1 = np.poly1d(z1)
p2 = np.poly1d(z2)
p3 = np.poly1d(z3)

# %%

print('p1 =\n', p1)
print('p2 =\n', p2)
print('p3 =\n', p3)

# %%


x1 = np.linspace(-2, 7, 100)
y1 = p1(x1)
y2 = p2(x1)
y3 = p3(x1)
plt.scatter(x, y)
plt.plot(x1, y1, label='linear')
plt.plot(x1, y2, label='quadratic')
plt.plot(x1, y3, label='cubic')
plt.legend()

输出如下：

p1 =
  
7.4 x - 4
p2 =
      2
0.5 x + 4.9 x - 1.5
p3 =
          3     2
-0.3333 x + 3 x - 0.6667 x + 2

有一次多项式 $7.4x-4$，二次多项式 $0.5x^2+4.9x-1.5$，三次多项式 $-0.3333x^3+3x^2-0.6667x+2$，曲线如图所示。