>>> from env_helper import info; info()
页面更新时间: 2023-06-24 12:30:55
运行环境:
    Linux发行版本: Debian GNU/Linux 12 (bookworm)
    操作系统内核: Linux-6.1.0-9-amd64-x86_64-with-glibc2.36
    Python版本: 3.11.2

2.1. SciPy-数值计算库

SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。 例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。 作为入门介绍,让我们看看如何用SciPy进行插值处理、信号滤波。

2.1.1. 最小二乘拟合

假设有一组实验数据(x[i], y[i]),我们知道它们之间的函数关系:y = f(x),通过这些已知信息,需要确定 函数中的一些参数项。例如,如果f是一个线型函数f(x) = k*x+b,那么参数k和b就是我们需要确定的 值。如果将这些参数用 p 表示的话,那么我们就是要找到一组 p 值使得如下公式中的S函数最小:

m ∑ S(p) = [y i − f (x i , p)] 2 i=1

这种算法被称之为最小二乘拟合(Least-square fitting)。 scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数leastsq。

下面是用leastsq进行数据拟合的一个例子:

>>> import numpy as np
>>> from scipy.optimize import leastsq
>>> import pylab as pl

定义数据拟合所用的函数: Asin(2pikx + theta) 。

>>> def func(x, p):
>>>     A, k, theta = p
>>>     return A*np.sin(2*np.pi*k*x+theta)

定义函数,实验数据x, y和拟合函数之间的差,p为拟合需要找到的系数:

>>> def residuals(p, y, x):
>>>     return y - func(x, p)
>>> x = np.linspace(0, -2*np.pi, 100)
>>> A, k, theta = 10, 0.34, np.pi/6 # 真实数据的函数参数
>>> y0 = func(x, [A, k, theta]) # 真实数据
>>> y1 = y0 + 2 * np.random.randn(len(x)) # 加入噪声之后的实验数据
>>>
>>>
>>> p0 = [7, 0.2, 0] # 第一次猜测的函数拟合参数

调用 leastsq 进行数据拟合, residuals 为计算误差的函数, p0 为拟合参数的初始值, args 为需要拟合的实验数据。

>>> plsq = leastsq(residuals, p0, args=(y1, x))
>>>
>>> print("真实参数:", [A, k, theta])
>>> print("拟合参数", plsq[0]) # 实验数据拟合后的参数
真实参数: [10, 0.34, 0.5235987755982988]
拟合参数 [10.74635075  0.34020524  0.46419834]

由于在制图中会用到中文,需要先设置字体:

>>> from matplotlib import font_manager
>>> font_manager.fontManager.addfont('/usr/share/fonts/truetype/wqy/wqy-zenhei.ttc')
>>> font_manager.fontManager.addfont('/usr/share/fonts/truetype/wqy/wqy-microhei.ttc')
>>> for font in font_manager.fontManager.ttflist:
>>>     # 查看字体名以及对应的字体文件名
>>>     if 'wqy' in font.name.lower() or 'wenquanyi' in font.name.lower():
>>>         print(font.name, '-', font.fname)
WenQuanYi Zen Hei - /usr/share/fonts/truetype/wqy/wqy-zenhei.ttc
WenQuanYi Micro Hei - /usr/share/fonts/truetype/wqy/wqy-microhei.ttc
>>> pl.rcParams['font.sans-serif'] = 'WenQuanYi Zen Hei'

进行可视化制图,调用leastsq函数对噪声正弦波数据进行曲线拟合:

>>> pl.plot(x, y0, label="真实数据")
>>> pl.plot(x, y1, label="带噪声的实验数据")
>>> pl.plot(x, func(x, plsq[0]), label="拟合数据")
>>> pl.legend()
>>> pl.show()
/usr/lib/python3/dist-packages/IPython/core/pylabtools.py:151: UserWarning: Glyph 8722 (N{MINUS SIGN}) missing from current font.
  fig.canvas.print_figure(bytes_io, **kw)
_images/sec01_square_17_1.png

这个例子中我们要拟合的函数是一个正弦波函数,它有三个参数 A, k, theta , 分别对应振幅、频率、相角。 假设我们的实验数据是一组包含噪声的数据 x, y1,其中y1是在真实数据y0的基础上加入噪声的到了。

通过leastsq函数对带噪声的实验数据x, y1进行数据拟合,可以找到x和真实数据y0之间的正弦关系的三个参数: A, k, theta。下面是程序的输出:

我们看到拟合参数虽然和真实参数完全不同,但是由于正弦函数具有周期性,实际上拟合参数得到的 函数和真实参数对应的函数是一致的。