数据处理基础——函数拟合与图像绘制

零 · 前言

笔者由于在刚刚过去的一个学期的大物实验报告中偷懒，图表全是用网格纸手绘的，现决定学习使用更现代化的数据处理与绘图方法以弥补当初错失的学习契机。此外，本文作为笔者的第一篇博客文章，也用于测试个性化前端页面下随笔内容的正常显示。

本文主要使用Python的Matplotlib、Numpy库，实现已知阶数的一元多项式函数拟合、已知函数结构的多元任意函数拟合、函数图像绘制。

一 · 已知阶数的一元多项式函数拟合

Numpy库的polyfit()函数基于最小二乘法对泰勒展开多项式拟合，给定多项式阶数，该函数可得出最佳拟合系数。

对于多项式函数\(f(x) = \sum_{i=0}^{n} a_{i} x^{i}\)，其误差平方和为\(E = \sum_{i=1}^{n}\left | f(x_{i}) - y_{i} \right |^{2}\)。

五阶多项式函数拟合举例

import numpy as np
import matplotlib.pyplot as plt

x = [30, 35, 40, 45, 50, 55, 60, 65, 70]
y = [2393, 1992, 1672, 1427, 1189, 1016, 892, 740, 632]

# 拟合n阶多项式
f = np.polyfit(x, y, 5)
print(np.poly1d(f)) # 由高次至低次输出多项式系数
y1 = [f[0] * i ** 5 + f[1] * i ** 4 + f[2] * i ** 3 + f[3] * i ** 2 + f[4] * i + f[5] for i in x]

# 绘图
plt.scatter(x, y, marker='.') # 散点图
plt.plot(x, y1) # 拟合曲线
plt.title('test')
plt.xlabel('para_x')
plt.xticks([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75])
plt.ylabel('para_y')
plt.yticks([0, 500, 1000, 1500, 2000, 2500])
plt.tight_layout()
plt.show()

二 · 已知函数结构的多元任意函数拟合

Scipy库的curve_fit()函数通过非线性最小二乘法拟合数据，给定函数结构，该函数可得出最佳拟合系数。

可用决定系数评估拟合结果\(R^2 = 1 - \frac{SSE}{SST}\)，其中残差平方和\(SSE = \sum_{i=1}^{n} (y_{i} - f(x_{i}))^{2}\)，离差平方和\(SST = \sum_{i=1}^{n} (y_{i} - \bar{y})^{2}\)。

二元函数拟合举例

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt

# 自定义函数结构
def func(para, a, b): 
    return a * para[0] ** 2 + b * para[1] ** 3

# 计算R²值
def estimate_r2(true, pred): 
    ss_tot = np.sum([(i - np.mean(true)) ** 2 for i in true])
    ss_res = np.sum([(i - j) ** 2 for i, j in zip(true, pred)])
    r2 = 1 - ss_res / ss_tot
    return r2

x = np.arange(0, 20)
y = np.arange(20, 40)
z = [i ** 2 + j ** 3 + np.random.uniform(-10, 10) for i, j in zip(x, y)]

# 拟合已知结构函数
p, cov = opt.curve_fit(func, [x, y], z) # p[i]为拟合系数，cov[i, j]为第i个系数与第j个系数的协方差
print(p) # 按序输出拟合系数
z_pred = [func([i, j], *p) for i, j in zip(x, y)]
r2 = estimate_r2(z, z_pred)
print(f"{r2:.10f}")

---

靡不有初，鲜克有终。浅开一坑，希望能坚持更新下去。