拟合曲线，其中每个点都是 ODE 的解

Question

我有以下代码：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import random

def dose(y, t, b, s, c, p, i):
    target, infectious, virus = y
    dydt = [-b*target*virus, b*target*virus - s*infectious, (1/(i+1))*p*infectious - c*virus]
    return dydt


b = 0.00001
s = 4
c = 4
p = 2000000
D = np.logspace(-3, 3, 7)
mylist = []

y0 = [1, 0, 0.01]
t = np.linspace(0, 60, 1000)
for i in D:
    sol = odeint(dose, y0, t, args=(b, s, c, p, i))
    #plt.plot(t, sol[:, 0], label='D = ' + str(i))
    V = sol[:, 2]
    mylist.append(V[48]/0.01950269536785707)


def add_noise(d, noise_pct):
    return [x + random.gauss(0, noise_pct * x) for x in d]

mat = np.array(mylist)

mylist2 = add_noise(mylist, 0.05)
mat2 = np.array(mylist2)
plt.plot(D, mat2)

plt.xscale('log')
#plt.plot(t, sol[:, 2], 'r', label='virus')
plt.legend(loc='best')
plt.xlabel('time')
plt.grid()
#plt.show()
print(D)
print(mat2)

现在 D 和 mat2 包含 x 和 y 值，我想用它们来拟合函数剂量中原始 ODE 系统的解。我也想绘制新的拟合函数。

我该怎么做？

为了澄清，我从原始图表中添加了 5% 到 7 个数据点的误差变化，并得到了一个 x 和 y 值的 np.array，我想找到一条最适合的曲线，它遵循与以下相同的微分方程原始图表。

Answer 1

让我们重新表述一下您的问题，使其更高效、更适合。

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize

我们编写您系统的动态以符合新的求解器接口

solve_ivp

:

def dose(t, y, b, s, c, p, d):
    target, infectious, virus = y
    return np.array([
        -b * target * virus,
        b * target * virus - s * infectious,
        (1. / (d + 1.)) * p * infectious - c * virus
    ])

现在我们编写的模型函数具有与

curve_fit

兼容的签名，我们还在求解ODE时进行了简化，以提高算法的时间复杂度：

def model(D, b, s, c, p):
    solutions = []
    for d in D:
        solution = integrate.solve_ivp(
            dose, [0, 5], y0=[1, 0, 0.01],
            t_eval=[2.8828828828828827], # Translating np.linspace(0, 60, 1000)[48]
            args=(b, s, c, p, d)
        )
        data = solution.y[2, 0] / 0.01950269536785707  # I'm intrigued by this part of the your code I translated
        solutions.append(data)
    return np.array(solutions)

现在我们生成一些合成数据集：

b = 0.00001
s = 4
c = 4
p = 2000000

D = np.logspace(-3, 3, 20)
z = model(D, b, s, c, p)
s = np.ones_like(z) * 0.01
n = s * np.random.normal(size=s.size)
zn = z + n

然后我们进行调整：

popt, pcov = optimize.curve_fit(
    model, D, zn, p0=[1e-5, 1, 1, 1e6],
    method="trf", bounds=(0, np.inf),
    sigma=s, absolute_sigma=True
)
#(array([1.88403946e-05, 3.44614223e+00, 4.26709794e+00, 1.00015953e+06]),
# array([[ 1.53054270e-12,  3.43500804e-05, -3.44589474e-05,
#         -1.09829485e-12],
#        [ 3.43500804e-05,  1.09331529e+03, -1.10009434e+03,
#         -3.49946049e-05],
#        [-3.44589474e-05, -1.10009434e+03,  1.10693945e+03,
#          3.52118585e-05],
#        [-1.09829485e-12, -3.49946049e-05,  3.52118585e-05,
#          1.12010298e-12]]))

参数数量级相似，适应度好：

Dlog = np.logspace(-3, 3, 200)
fig, axe = plt.subplots()
axe.scatter(D, zn)
axe.semilogx(Dlog, model(Dlog, *popt))
axe.grid()

但是这种拟合对噪声和初始猜测非常敏感，这是预期的，因为它必须为曲线的每个点求解刚性 ODE。参数的多种组合似乎也可以导致等效的适应度。