原文：
towardsdatascience.com/differentiate-noisy-time-series-data-with-symbolic-regression-2e27db2392c8注意：如果你没有 Medium 订阅，你可以在这里免费阅读文章 here！
时间序列轮廓在我们的日常生活中无处不在。也有许多专门的研究工作处理这些问题。
简而言之，时间序列轮廓是一系列后续数据点 y(0), y(1), … ,y(t) 的集合，其中在时间 t 的一个点依赖于时间 t-1 的前一个点（甚至更早的时间）。
在许多应用中，人们感兴趣的是在有一些先前点的情况下预测轮廓的行为。为此，存在各种各样的建模方法。在核心上，这些模型可能需要一些关于过去（或现在）的信息，并给出关于未来轮廓外观的估计。人们可以找到很多处理此类时间序列预测的工作，例如使用神经网络描述天气（Bi et al., 2023），通过深度学习分析股价行为（Xiao and Su, 2022），或者分析药品需求演变（Rathipriya et al., 2023)。当然，这些研究工作是我快速搜索后找到的，所以还有很多其他的东西存在。
然而，正如你所知，为了使模型能够预测时间序列轮廓，我们需要数据。通常，数据越好，我们就能更好地描述研究的过程。
如果我们给它一些旧的时间点 y(0), … , y(t)，人们可以训练一个模型来预测下一个状态 y(t+1)。然而，在某些应用中，我们可能需要一个模型，它接受当前的观察结果 y(0), … , y(t)，并预测系统在当前或下一个时间点的变化速度。因此，我们想要系统的导数 dy 而不是可观察的状态 y。因此，为了训练这样的模型并使其输出导数 dy，我们首先需要收集这些导数。通常，这些导数是直接从观察到的数据 y 中计算出来的，因为根据情况，测量导数可能很困难，甚至不可能。
这就是一个小问题让我们的生活变得如此困难的地方：噪声。
每个处理过噪声时间序列数据的人都知道计算它们的导数可能有多么痛苦。
处理时间序列数据中的噪声有许多可能性。在这篇文章中，我想向您展示一种在我们的一个研究项目中工作得相当好的可能性，在该项目中，时间序列轮廓中可用的数据点并不多。
所以，给自己倒一杯咖啡 ☕，然后启动 Python 🐍！
数据生成
要推导数据，我们需要数据。让我们创建一个场景。正如提到的，我想创建一个我们没有数千个数据点的场景。
为了做到这一点，让我们考虑一个反应堆，其中一些生物质（即，细菌或其他细胞）消耗一个底物（即，一些糖——对我们来说，它就像咖啡一样，对它们的生命至关重要 ☕），并产生一个产品，例如一种作为药物有用的蛋白质。因此，我们有“生物质”，“底物”，以及我们的“产品”。我们可以在特定时间点 t 直接以每升克（g/L）为单位测量这些物种。所以它们在反应器中的浓度基本上是我们的可观察数据 y。
我上面提到的研究工作的总体目标是创建一个特定形式的模型，该模型输出导数。所以这里我不深入细节，但这篇文章的重点是我们需要推导出观察到的时序轮廓 y。
所以让我们用 Python 来做这件事，首先导入一些包。此外，我们定义了一个用于后续使用的图形保存函数。
import os
import random
random.seed(0)
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from sklearn.metrics import mean_squared_error
# Function to save figures
FIG_SIZE = (7, 2.5)
FIG_DPI = 400
FIG_SAVEPATH = '.'
FIG_SAVEEXT = ['png']
def save_figure(figure_savepath : str = None,
save_figures_extension : list = ['png'],
figure_name : str = 'Figure',
dpi : int = 1200):
figure_name = figure_name
if isinstance(save_figures_extension, list):
figure_extension_list = save_figures_extension
else:
raise ValueError('[-] The indicated file extension for figures needs to be a list!')
for figure_extension in figure_extension_list:
if figure_savepath is None:
figure_savepath = os.getcwd()
savepath = os.path.join(figure_savepath, figure_name+'.'+figure_extension)
plt.savefig(savepath, dpi=dpi)
print(f'[!] Figure saved as: {savepath}')
接下来，我们可以定义一个类来存储我提到的场景。在这里，我不会深入生物系统的细节，但如果你有任何问题，请随时在评论中提问 😊 代码将使用给定的初始条件 y0 和 80 小时的时间跨度（tspan）以及 20 个时间点来求解常微分方程组（ODEs）。正如你所看到的，采样频率因此不是很高，大约每四小时一个样本。
注意：常微分方程组（ODE）系统是从我在函数描述中提到的来源改编的。这个例子（以及其他例子）在 Python 中是现成的，可以通过insidapy（也可以通过 pip 安装）获得。
#%% Generate time series data
# ============================================================================
class datagen():
# ------------------------------------------------------------------------
def __init__(self,
y0 = np.array([0.1, 60, 0]),
tspan = np.linspace(0, 80, 20)):
self.y0 = y0
self.tspan = tspan
self.specnames = ['Biomass', 'Substrate', 'Product']
# ------------------------------------------------------------------------
def solve_ODE_model(self):
'''Solve the ODE model for the given sample and time span.
Returns the derivatives (dydt) as an array of shape [n,] and the runtime of the ODE solver.
'''
# Simulate the reactor operation until the selected time tf
self.y = odeint(func=self.ODEmodel, y0=self.y0, t=self.tspan)
# ------------------------------------------------------------------------
def ODEmodel(self, y, t):
'''Adjusted ODE model for batch fermentation. Literature: Turton, Shaeiwtz, Bhattacharyya, Whiting
"Analysis, synthesis and design of chemical processes", Prentice Hall. 2018.
'''
# Variables
X = y[0]
S = y[1]
P = y[2]
# Parameters
mu_max = 0.25; #h^-1
K_S = 105.4; #kg/m³
Y_XS = 0.07; #[-]
Y_PS = 0.167; #[-]
KXmu = 121.8669;#g/L constant for inhibition of biomass growth caused by higher cell densities
T = 273 + 35; #K
R = 0.0083145; #kJ/(K*mol) universial gas constant
k1_ = 130.0307; #[-] constant for activation of biomass growth
E1_ = 12.4321; #kJ/mol activation enthalpy for biomass growth
k2_ = 3.8343e48; #[-] constant for inactivation of biomass growth
E2_ = 298.5476; #kJ/mol inactivation enthalpy for biomass growth
# Define temperature dependency of the rate constants
k1 = k1_ * np.exp(-E1_ /(R*T))
k2 = k2_ * np.exp(-E2_ /(R*T))
# Calculate specific growth rate of the biomass
mu = (mu_max*S)/(K_S+S) * k1/(1+k2) * (1-(X/(KXmu+X)))
# Calculate consumption of substrate
sigma = -(1/Y_XS)*mu
# Calculate specific production rate of the protein
pi = Y_PS/Y_XS*mu
# Vectorization of rates
rate = np.hstack((mu.reshape(-1,1), sigma.reshape(-1,1), pi.reshape(-1,1)))
# ODEs of biomass, volume and product
dydt = rate * X.reshape(-1,1)
# Return
return dydt.reshape(-1,)
# ------------------------------------------------------------------------
def addnoise_per_species(self, percentage=5):
'''Uses some ground truth data and adds some noise to it.
'''
self.y_noisy = np.zeros((self.y.shape[0], 0))
for spec_id in range(self.y.shape[1]):
rmse = mean_squared_error(self.y[:, spec_id], np.zeros(self.y[:, spec_id].shape), squared=False)
y_noisy_spec = self.y[:, spec_id] + np.random.normal(0, rmse / 100.0 * percentage, self.y[:, spec_id].shape)
self.y_noisy = np.hstack((self.y_noisy, y_noisy_spec.reshape(-1,1)))
# ----------------------------------------
def evaluate_true_derivatives(self):
'''Evaluate the true derivatives of the generated data.
'''
self.y_true_diff = np.zeros((0, self.y.shape[1]))
for t_id, t in enumerate(self.tspan):
self.y_true_diff = np.vstack((self.y_true_diff, self.ODEmodel(y=self.y[t_id, :], t=t).reshape(1,-1)))
有了这个类，我们可以实例化它，通过求解常微分方程组（ODE）来创建数据，并使用噪声破坏获得轮廓。
很酷，所以查看图 1，你可以看到真实数据（即，系统在没有噪声的情况下会是什么样子）作为虚线黑线，以及噪声数据作为黑点。所以你看，生物质正在增长（左图——细菌吃并生长），糖被消耗（中图——细菌吃它），产品正在形成（右图——它们产生蛋白质）：
# Instantiate the class
data = datagen()
# Create data and add noise
data.solve_ODE_model()
data.addnoise_per_species()
data.evaluate_true_derivatives()
# Plot true and noisy data
fig, ax = plt.subplots(figsize=FIG_SIZE, ncols=3)
for i in range(len(data.specnames)):
ax[i].plot(data.tspan, data.y[:, i], marker='', linestyle='--', color='black', label='Ground truth')
ax[i].plot(data.tspan, data.y_noisy[:, i], marker='o', markersize=4, linestyle='', color='black', alpha=0.6, label='Observed')
ax[i].set_xlabel('Time / h', fontsize=8)
ax[i].set_ylabel('{} / g/L'.format(data.specnames[i]), fontsize=8)
ax[i].tick_params(axis='both', which='major', labelsize=8)
if i == 1:
ax[i].legend(frameon=False, fontsize=8, loc='lower left')
plt.tight_layout()
save_figure(figure_savepath = FIG_SAVEPATH, save_figures_extension = FIG_SAVEEXT, figure_name = 'True_and_noisy_data', dpi = FIG_DPI)
图 1. 研究系统的基础真实情况（虚线）以及观察到的噪声数据点（黑色圆圈）。图由作者使用上述代码创建。图由作者创建。
由于我们想要计算这些轮廓的导数，因此拥有“真实”导数会很好，这样就可以比较所考虑的方法的表现。实际上，我们这里也有这些“真实”导数，它们基本上是上述代码中找到的常微分方程（ODE）系统的右侧。
因此，有了这些，我们有了我们的观测数据。现在让我们比较一些导数计算方法，看看它们的性能如何。
直接应用有限差分
在第一部分，让我们直接将有限差分（FD）方法应用于噪声数据。这是一种计算函数导数的简单方法。它定义如下：
其中 f’(x) 是函数 f(x) 在点 x 处的导数，而 h 是一个很小的步长。
让我们用我们的数据来做这件事（不是最有效的方法，但我希望代码是易于理解的）：
# Function to perform the finite difference given some x and y data
def finite_difference(x_data, y_data):
# Create empty vector for the finite difference approach (x and y vectors)
x_diff = np.zeros((y_data.shape[0]-1, y_data.shape[1]))
y_diff= np.zeros((y_data.shape[0]-1, y_data.shape[1]))
# Calculate finite difference derivatives
dt = (x_data[1] - x_data[0])/2
for spec_id in range(y_data.shape[1]):
# Initialize empty lists for derivatives
y_diff_FD = []
x_diff_FD = []
for t_ in range(len(x_data)-1):
x_diff_FD.append(x_data[t_] + dt)
y_diff_FD.append((y_data[t_+1, spec_id] - y_data[t_, spec_id])/(x_data[t_+1] - x_data[t_]))
# Convert lists to arrays
x_diff[:, spec_id] = np.array(x_diff_FD)
y_diff[:, spec_id] = np.array(y_diff_FD)
return x_diff[:,0].reshape(-1,), y_diff
# Perform finite difference
data.x_diff_FD, data.y_diff_FD = finite_difference(data.tspan, data.y_noisy)
那么，就是这样，让我们看看这些导数看起来如何，并将 FD 导数与“真实”导数进行比较：
fig, ax = plt.subplots(figsize=FIG_SIZE, ncols=3)
for i in range(len(data.specnames)):
ax[i].plot(data.tspan, data.y_true_diff[:, i], marker='', linestyle='--', color='black', label='True')
ax[i].plot(data.x_diff_FD, data.y_diff_FD[:, i], marker='', markersize=4, linestyle='-', color='blue', label='FD')
ax[i].set_xlabel('Time / h', fontsize=8)
ax[i].set_ylabel('Derivative {} / g/L/h'.format(data.specnames[i]), fontsize=8)
ax[i].tick_params(axis='both', which='major', labelsize=8)
if i == 1:
ax[i].legend(frameon=False, fontsize=8, loc='upper center')
plt.tight_layout()
save_figure(figure_savepath = FIG_SAVEPATH, save_figures_extension = FIG_SAVEEXT, figure_name = 'Derivative_comparison_FD', dpi = FIG_DPI)
图 2. 基础真实导数（虚线）与通过 FD 计算得到的导数（蓝色线条）。图由作者使用上述代码创建。图由作者创建。
在图 2 中，真实导数以虚线表示，而使用 FD 方法计算的导数以蓝色实线表示。正如你所见，有限差分方法对噪声非常敏感。因此，在这里应用有限差分方法并不是很有用，因为它会导致噪声放大，这意味着我们无法使用得到的数据进行任何操作。
过滤技术：Savitzky-Golay
处理噪声的一种方法是在我们处理数据之前尝试将其过滤掉。Savitzky-Golay（SG）滤波器是一种非常著名的平滑数据的方法。当然，还有许多其他可用的方法，但让我们考虑这种方法作为此类过滤技术的代表性方法。SG 滤波器在移动窗口中对数据进行多项式拟合。它在scipy包中很容易获得，因此让我们使用给定的窗口大小和多项式阶数将此滤波器应用于我们的数据：
注意：当然，窗口大小和阶数是可调整的超参数。但在此阶段，我们保持它们如下所述。
# Import SG filter
from scipy.signal import savgol_filter
# Define SG filter settings
window_size = 5
polynomial_order = 2
# Smooth the noisy data for each species
y_smooth = np.zeros((data.y_noisy.shape[0], data.y_noisy.shape[1]))
for i in range(data.y_noisy.shape[1]):
y_smooth[:, i] = savgol_filter(data.y_noisy[:, i], window_size, polynomial_order)
# Plot the smoothed data together with the noisy data
fig, ax = plt.subplots(figsize=FIG_SIZE, ncols=3)
for i in range(len(data.specnames)):
ax[i].plot(data.tspan, data.y[:, i], marker='', linestyle='--', color='black', label='True')
ax[i].plot(data.tspan, data.y_noisy[:, i], marker='o', markersize=4, linestyle='', color='black', label='Observed')
ax[i].plot(data.tspan, y_smooth[:, i], marker='', linestyle='-', color='red', label='Smoothed')
ax[i].set_xlabel('Time / h', fontsize=8)
ax[i].set_ylabel('{} / g/L'.format(data.specnames[i]), fontsize=8)
ax[i].tick_params(axis='both', which='major', labelsize=8)
if i == 1:
ax[i].legend(frameon=False, fontsize=8, loc='lower left')
plt.tight_layout()
save_figure(figure_savepath = FIG_SAVEPATH, save_figures_extension = FIG_SAVEEXT, figure_name = 'Smoothing_SG', dpi = FIG_DPI)
图 3. 基准真实数据（虚线黑色线）以及观测到的噪声数据点（黑色圆圈），以及 SG 平滑后的轮廓（实线红色线）。图由作者使用上述代码创建。图由作者创建。
太好了，让我们使用平滑后的数据点再次使用 FD 方法计算导数（基本上使用上面的相同函数）：
# Savitzky-Golay-smoothed data for differentiation with FD
data.x_diff_SGFD, data.y_diff_SGFD = finite_difference(data.tspan, y_smooth)
我们还可以将这部分可视化，以及“真实”的导数：
# Visualize true derivatives, FD and SGFD for each species
fig, ax = plt.subplots(figsize=FIG_SIZE, ncols=3)
for i in range(len(data.specnames)):
ax[i].plot(data.tspan, data.y_true_diff[:, i], 'k--', label='True')
ax[i].plot(data.x_diff_FD, data.y_diff_FD[:, i], 'b-', label='FD')
ax[i].plot(data.x_diff_SGFD, data.y_diff_SGFD[:, i], 'r-', label='SGFD')
ax[i].set_xlabel('Time / h', fontsize=8)
ax[i].set_ylabel('Derivative {} / g/L/h'.format(data.specnames[i]), fontsize=8)
ax[i].tick_params(axis='both', which='major', labelsize=8)
if i == 1:
ax[i].legend(frameon=False, fontsize=8, loc='best')
plt.tight_layout()
save_figure(figure_savepath = FIG_SAVEPATH, save_figures_extension = FIG_SAVEEXT, figure_name = 'Derivative_comparison_FD_SGFD', dpi = FIG_DPI)
图 4. 基准真实导数（虚线黑色线）以及通过 FD 计算得到的导数（实线蓝色线），以及使用 FD 平滑和导数的数据（实线红色线）。图由作者使用上述代码创建。图由作者创建。
因此，我们可以看到（使用 SG 滤波器的给定参数），我们似乎可以在导数计算上得到一些积极的效果。然而，这种方法并不完美，导数仍然相当噪声。让我们看看一种——根据我的经验——不太常用的方法。
另一个选项：符号回归
如果我们有大量的数据点，SG 滤波器可能会给出非常好的结果，因为平滑非常有效。然而，在我们的例子中，这似乎很难。
另一种处理噪声的方法是首先拟合一个函数到数据上。这个函数需要是我们所知的，这样我们就可以解析地推导这个函数。这种方法的简单方式是通过数据拟合一个多项式。这给我们一个以下形式的拟合函数：
其中 a0, a1, … , an 是多项式的系数，这些系数是通过拟合过程得到的。现在我们有一个封闭形式的函数可用，我们可以对这个函数进行解析推导。有了这个，我们也得到了导数的解析表达式，我们可以在给定的时间点进行评估。
这种方法的缺点是可能很难找到一个好的多项式阶数，或者一般的功能形式。如果阶数太低，多项式可能无法很好地拟合数据。如果阶数太高，多项式可能拟合噪声。
如果我们有数据可用，有一个工具可以为我们创建一个封闭形式的表达式：符号回归（SR）。简而言之，SR 寻找可以很好地解释我们数据的方程。更合适的方程被存储起来，而那些不能很好地代表我们数据的方程则被忽略。我们只需要告诉算法（除了我们提供的数据外），最终方程中允许哪些运算符（即乘法、加法、除法、指数等）。
我在这里不深入介绍该方法，但如果您对这个方法的工作原理或应用感兴趣，您可以在 Medium 上查看我关于这个主题的其他文章：
您还可以查看可用的代码文档（例如 PySR——一个超级棒的工具！），或者一些已发表的科研工作（yebbo，其中之一是我们的，其他的是比我聪明得多的人 😁）：Guimerà 等人 (2020)，Forster 等人 (2024)，de Carvalho Servia (2024)。
注意：接下来的部分将涉及许多步骤，以及相当多的代码。正如之前提到的，我们在我们的一篇研究文章中应用了这种方法，您可以在这里找到它（开放获取）。
它基于由 Guimerà 等人 (2020) 创建的免费可用的符号回归算法——他们将其命名为贝叶斯机器科学家（BMS）。我们考虑到微分任务的目的，将该算法打包成 Python 包 udiff，并在 GitHub 上也使其可用。
我们在代码中包含了一些停止条件（您也可以在下面的代码片段中看到）：
-
nsteps：算法应运行的步数。
-
maxtime：算法应运行的最大时间（以秒为单位）。
-
minr2：算法应达到的最小 R2 值（确定系数）。
然后，休息一下，再来点咖啡，这对我们的生活来说同样至关重要 ☕
准备好了吗？😊 所以，我们有一些噪声数据，我们希望符号回归算法给我们一个方程，该方程可以“平滑”数据（使用 udiff.smooth）。然后，我们希望对这个方程进行解析微分（使用 udiff.differentiate）。让我们来做这件事：
注意：这个过程可能需要很长时间，所以请坐下来享受您的咖啡 ☕ 它也是一个随机过程。所以每次执行拟合时，它可能会导致不同的结果！
# Import codes
from udiff.smooth import smooth_bms # Smooths the data using a symoblic regression model
from udiff.differentiate import differentiator # Differentiates the obtained symbolic regression model analytically
# Create empty arrray for fitted profiles
data.y_smooth_sr = np.zeros(data.y.shape)
data.y_diff_sr = np.zeros(data.y.shape)
# Fit profile and derive it for each species (takes some time!)
for spec_id in range(data.y.shape[1]):
# Each species individually
X = data.tspan
Y = data.y_noisy[:, spec_id]
# Create new object
obj = smooth_bms(x=X, y=Y, scaling=False)
# Fit profiles
obj.fit_bms(nsteps=1e4, maxtime=1800, minr2=0.999, show_update=True, update_every_n_seconds=200)
# Store data
data.y_smooth_sr[:, spec_id] = obj.y_smooth
# Differentiate algebraic equation
diffobj = differentiator(obj)
diffobj.differentiate()
# Store data
data.y_diff_sr[:, spec_id] = diffobj.y_diff
在此过程中，您应该会得到一些输出。
在幕后发生的事情是，SR 工具正在寻找一个很好地拟合我们数据的函数。如果我们绘制它，我们会得到与图 3 中所示相似的图像。美丽之处在于我们确切地知道这个函数的样子，因此我们可以对其进行解析。udiff包会自动为我们完成这项工作，依赖于sympy的功能。
因此，在完成这些之后，让我们可视化 SR 导数以及之前通过 FD 和 SG-FD 方法计算得到的导数：
# Visualize true derivatives, FD, SGFD and SR for each species
fig, ax = plt.subplots(figsize=FIG_SIZE, ncols=3)
nm = data.specnames
for i in range(len(data.specnames)):
ax[i].plot(data.tspan, data.y_true_diff[:, i], 'k--', label='True')
ax[i].plot(data.x_diff_FD, data.y_diff_FD[:, i], 'b-', label='FD')
ax[i].plot(data.x_diff_SGFD, data.y_diff_SGFD[:, i], 'r-', label='SGFD')
ax[i].plot(data.tspan, data.y_diff_sr[:, i], 'g-', label='SR')
ax[i].set_xlabel('Time / h', fontsize=8)
ax[i].set_ylabel('Derivative {} / g/L/h'.format(data.specnames[i]), fontsize=8)
ax[i].tick_params(axis='both', which='major', labelsize=8)
if i == 1:
ax[i].legend(frameon=False, fontsize=8, loc='best')
plt.tight_layout()
save_figure(figure_savepath = FIG_SAVEPATH, save_figures_extension = FIG_SAVEEXT, figure_name = 'Derivative_comparison_FD_SGFD_BMS', dpi = FIG_DPI)
图 5. 基础真实导数（虚线黑色线条）以及通过 FD 计算得到的导数（实线蓝色线条）、使用 FD 平滑和导出的数据（实线红色线条），以及符号回归平滑后的数据随后进行解析微分（实线绿色线条）。图由作者使用上述代码创建。图由作者创建。
如您所见，SR 方法给出了相当不错的结果。导数相当平滑，并且对导数计算的噪声影响得到了降低。现在，由于这些数据点与“真实”导数更接近，基于这些数据点训练的模型很可能也会导致更好的预测。
当然，像其他方法一样，这个方法也不是完美的。然而，它可能是一个很好的、有用的方法，可以添加到您的工具箱中 😊
结论
我希望我能向您展示拥有封闭形式方程的强大之处，以便找到导数的良好近似。当然，这并不是唯一的方法，但我认为它确实是一个替代方案，并且在特定情况下非常有用。
对一些人来说可能有所帮助（如果确实如此，请查看我们的研究工作[https://www.sciencedirect.com/science/article/pii/S0009250924009060]，我们在那里更详细地讨论了所有这些）。即使它对您的工作没有直接的帮助，我也希望您能享受阅读这个故事，并且在阅读的同时至少能喝上一杯美味的咖啡 😊 ☕️
就这样，感谢阅读并祝大家愉快地学习！😊
对于那些直接跳到结论的人来说，这里还有一个故事，我将在其中讨论符号回归与优化。如果您感兴趣，请查看它——当然——如果您还有足够的咖啡继续在 Medium 上阅读 ☕️
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
