于散乱中发掘纯真——拟合算法

---

联系数计加油站获取论文和代码等资源。
邮箱:MathComputerSharing@outlook.com
微信公众号:数计加油站
致谢:清风数学建模

---

规律是纯真的，数据是散乱的😊

我们的世界所遵从的规律，往往是简单而纯真的。我们可以使用准确的数学语言描述出一条条优美的规律。比如物理世界中的杠杆原理，万有引力定律，质能方程，相对论…

我们深信，万事万物都有它背后隐藏的规律。而且这种规律，往往比较直白简单。

发掘这些规律的第一步，是通过实验或者说实践，获取到一系列相关数据。有了数据，我们就可以利用以往的经验，对数据背后的规律进行一定的猜测。

不过，有个小小的问题：数据在获取的过程中容易产生误差，而且这种误差是永远不可能避免，只能尽量减小的。从包含误差的数据中挖掘简单的规律，这是技术含量很高的一件事。

所以说，我们猜测规律，往往只需要追求规律与数据的神似，而不需要规律能生成每个数据。

有了规律，开始拟合😁

猜测规律，就是写出一个带参数的函数表达式。比如二维数据，我们猜测它符合直线规律，我们就会猜测数据接近于一条直线 $y=kx+b$；我们猜测它符合抛物线规律，我们就会猜测数据接近于一条抛物线$y=a{x}^2+bx+c$;我们猜测它符合正弦规律，我们就会猜测数据接近于一条正弦曲线$y=Asin(\omega x+\phi )$

怎么得到确切的参数值呢?这个过程就称为拟合。

拟合要求数据离规律最接近。这就需要引入数据点与拟合数据点的距离的概念。常用的距离衡量式是最小二乘法中使用到的
$$d_i=({\hat{y}}_i-y_i{)}^2$$

拟合的目的就是在已划定范围的函数族中，找到能够使所有数据点与拟合数据点距离之和最小的参数。有了这些参数，我们猜测的那个规律就完全确定下来了。

哪个规律猜的好😎

可能出现这样的情况。张三猜测的规律是这样的，李四猜测的规律却是那样的，看着都很有道理的样子。怎样来评价谁的拟合更好呢？

评价准则一般有两个：
（1）接近程度。即已知数据点与拟合数据点越接近越好。相关的评价指标有$SSE$和$R^2$等
（2）复杂度。即接近程度差不多的时候，越简单的规律表达式往往越好，越get到规律的本质

一篇简短的小论文

Matlab 代码

自定义的线性拟合函数

function [theta,SSE,R_square] = LinearRegression(x,y,functions)
% 线性拟合
% x,y均为列向量
% functions为元胞函数数组
    
    %% 构造系数矩阵
    m = length(x);
    n = length(functions);
    R = ones(m,n);
    for i = 1:n
        func  = functions{i};
        R(:,i) = func(x);
    end
    
    %% 解正规方程
    A = (R' * R);
    b = R' * y;
    theta =  A \ b;
    
    %% 计算y的拟合值
    hat_y = ones(m,1);
    for i = 1:m
        gi = 0;
        for j = 1:n
        gi = gi+theta(j)*functions{j}(x(i));
        end
        hat_y(i) = gi;
    end
    
    %% 计算SSE和R-square
    mean_y = mean(y);
    SSE = sum((hat_y-y).^2);
    SSR = sum((hat_y-mean_y).^2);
    SST = sum((y-mean_y).^2);
    R_square = SSR/SST;
end

一次函数和二次函数拟合比较

clear,clc

%% 输入数据
x = 0:7;
y = [27;26.8;26.5;26.3;26.1;25.7;25.3;24.8];

%% 输入拟合函数
f0 = @(x) 1;
f1 = @(x) x;
f2 = @(x) x.^2;
functions1 = {f0,f1};
functions2 = {f0,f1,f2};

%% kx+b拟合
[theta1,SSE1,R_square1] = LinearRegression(x',y,functions1);

%% ax^2+bx+c拟合
[theta2,SSE2,R_square2] = LinearRegression(x',y,functions2);

%% 可视化
new_x = linspace(min(x),max(x),1000);
new_y1 = ones(size(new_x));
new_y2 = ones(size(new_x));
for i = 1:length(new_x)
    gi1 = 0;
    for j = 1:length(functions1)
        gi1 = gi1+theta1(j)*functions1{j}(new_x(i));
    end
    new_y1(i) = gi1;
    gi2=0;
    for j =1:length(functions2)
        gi2 = gi2+theta2(j)*functions2{j}(new_x(i));
    end
    new_y2(i) = gi2;
end

plot(x,y,'o',new_x,new_y1,new_x,new_y2)
legend("原始数据点","kx+b","ax^2+bx+c")

Python代码

初始化

# 导包
import numpy as np
import matplotlib.pyplot as plt

np.set_printoptions(suppress=True)

plt.rcParams["font.sans-serif"] = ["SimHei"]  #设置字体
plt.rcParams["axes.unicode_minus"] = False  # 解决图像中的“-”负号的乱码问题

# 初始化数据
x = np.arange(0, 8)
y = np.array([27, 26.8, 26.5, 26.3, 26.1, 25.7, 25.3, 24.8])

# 初始化函数
f0 = lambda x: 1
f1 = lambda x: x
f2 = lambda x: x ** 2

functions1 = [f0,f1]
functions2 = [f0,f1,f2]

函数定义

def LinearRegression(x, y, functions):
"""线性拟合"""

    # 参数初始化
    m = len(x)
    n = len(functions)
    R = np.ones((m, n))

    for i in range(n):
        func = functions[i]
        R[:, i] = func(x)

    # 解正规方程
    A = np.dot(R.T, R)
    b = np.dot(R.T, y)
    theta = np.dot(np.linalg.inv(A), b)

    # 计算拟合函数值
    hat_y = np.ones(m)
    for i in range(m):
        gi = 0
        for j in range(n):
            gi = gi + theta[j] * functions[j](x[i])
        hat_y[i] = gi

    # 计算SSE和R_square
    mean_y = np.mean(y)
    SSE = np.sum((hat_y - y) ** 2)
    SSR = np.sum((hat_y - mean_y) ** 2)
    SST = np.sum((y - mean_y) ** 2)
    R_square = SSR / SST

    return theta, SSE, R_square

一次函数与二次函数拟合对比

# kx+b拟合
theta1,SSE1,R_square1 = LinearRegression(x,y,[f0,f1])

# ax^2+bx+c拟合
theta2,SSE2,R_square2 = LinearRegression(x,y,[f0,f1,f2])

# 可视化
new_x = np.linspace(np.min(x),np.max(x),1000)
new_y1 = np.ones(new_x.shape)
new_y2 = np.ones(new_x.shape)
for i in range(len(new_x)):
    gi1 = 0
    for j in range(len(functions1)):
        gi1 = gi1 + np.dot(theta1[j],functions1[j](new_x[i]))
    new_y1[i] = gi1

    gi2 = 0
    for j in range(len(functions2)):
        gi2 = gi2 + np.dot(theta2[j],functions2[j](new_x[i]))
    new_y2[i] = gi2

plt.plot(x,y,'o',new_x,new_y1,new_x,new_y2)
plt.legend(["原始数据点","$kx+b$","$ax^2+bx+c$"])