PyMC与Sunode中时变ODE的符号插值实现问题求助

针对你在sunode/PyMC场景中，符号类型t无法实现输入数据插值的问题，以下是几种可行的解决思路：

1. 用SymPy拟合符号插值函数 #

将离散的输入数据拟合为可用于符号t的插值函数，比如样条插值。SymPy的interpolate函数可直接生成符号表达式，能被sunode等求解器处理：

import sympy as sp
import numpy as np

# 假设influent包含时间点t_in和对应值y_in
t_in = np.array([0, 1, 2, 3])
y_in = np.array([5, 7, 6, 8])

# 创建符号变量
t_sym = sp.Symbol('t')

# 拟合三次样条的符号表达式
influent_sym_spline = sp.interpolate(list(zip(t_in, y_in)), t_sym)

# 在ODE函数中使用该符号表达式
def ode(y, t, theta, influent_sym):
    # 代入符号t计算输入值
    inflow = influent_sym.subs('t', t)
    # 后续ODE逻辑示例
    dy_dt = theta * y + inflow
    return dy_dt

2. 构建分段符号插值函数 #

若样条插值复杂度较高，可手动构建分段线性的符号函数，通过sympy.Piecewise实现轻量化的符号插值：

def create_piecewise_linear(t_in, y_in):
    t_sym = sp.Symbol('t')
    pieces = []
    # 生成每个区间的线性插值表达式
    for i in range(len(t_in)-1):
        t_start, t_end = t_in[i], t_in[i+1]
        y_start, y_end = y_in[i], y_in[i+1]
        slope = (y_end - y_start)/(t_end - t_start)
        expr = y_start + slope*(t_sym - t_start)
        pieces.append((expr, sp.And(t_sym >= t_start, t_sym < t_end)))
    # 处理区间外的边界情况
    pieces.append((y_in[-1], t_sym >= t_in[-1]))
    pieces.append((y_in[0], t_sym < t_in[0]))
    return sp.Piecewise(*pieces)

# 生成分段函数
influent_piecewise = create_piecewise_linear(t_in, y_in)

# 在ODE中使用
def ode(y, t, theta, influent_piecewise):
    inflow = influent_piecewise.subs('t', t)
    # ODE逻辑示例
    dy_dt = theta * y + inflow
    return dy_dt

3. 拆分积分区间适配求解器 #

若符号插值性能不足，可利用Sundials求解器的分段积分特性，避开符号插值问题：

• 将观测时间序列拆分为输入数据点之间的子区间
• 对每个子区间，用数值插值获取输入值并求解ODE
• 拼接各子区间的结果保证状态连续性

在PyMC中配合sunode的实现示例：

import pymc as pm
import sunode
import numpy as np

# 观测时间点与数据
obs_times = np.linspace(0, 3, 10)
obs_data = np.random.normal(size=10)
# 输入数据的时间点
t_in = np.array([0, 1, 2, 3])
y_in = np.array([5, 7, 6, 8])

with pm.Model() as model:
    theta = pm.Normal('theta', 0, 1)
    y0 = pm.Normal('y0', 0, 1)
    
    y_prev = y0
    y_vals = []
    # 遍历每个输入区间求解ODE
    for i in range(len(t_in)-1):
        t_start, t_end = t_in[i], t_in[i+1]
        # 筛选当前区间内的观测点
        interval_times = obs_times[(obs_times >= t_start) & (obs_times <= t_end)]
        if len(interval_times) == 0:
            continue
        # 取区间中点的插值输入值
        inflow_val = np.interp((t_start + t_end)/2, t_in, y_in)
        # 求解当前区间的ODE
        sol = sunode.solve_ivp(
            y0={'y': y_prev},
            params={'theta': theta},
            times=interval_times,
            rhs=lambda y, t, p: {'y': p['theta'] * y['y'] + inflow_val},
        )
        y_vals.append(sol['y'])
        y_prev = sol['y'][-1]
    
    # 拼接结果并拟合观测数据
    y_obs = pm.Normal('y_obs', mu=np.concatenate(y_vals), sigma=0.1, observed=obs_data)
    idata = pm.sample()

关键注意事项 #

• 符号插值函数需保证可微分，确保PyMC采样能计算梯度（SymPy样条、分段线性函数均满足要求）
• 使用sunode时，避免在符号表达式中引入SymPy不支持转换为C代码的函数
• 分段积分需保证子区间边界的状态连续性，避免求解结果出现断层

内容的提问来源于stack exchange，提问作者Rvdb