在给定其他条件的情况下，在特定时间点估计 ODE 中的参数

Question

假设我的 ODE 系统中只有一个参数。我想推断这一点。我是否必须简单地重新排列等式以隔离值？在有多个方程式的系统中，这是如何完成的？例如：

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

#three compartments, Susceptible S, infected I, recovered R
#dS/dt, dI/dt, dR/dt
#susceptible sees birth rate coming in, deaths leaving and force of infection leaving
#infected sees FOI coming in, deaths leaving and recovery rates
#recovered sees recovery rate coming in, deaths leaving
#beta is tranmission coefficient, FOI is beta * (I/N) where N is total pop
#initially consider a model not accounting for births and deaths




# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0 = 10, 0
# Everyone else, S0, is susceptible to infection initially.
S0 = N - I0 - R0
J0 = I0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
#beta = 
gamma = 1/7
# A grid of time points (in days)
t = np.linspace(0, 100, 100+1)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma):
    S, I, R, J = y
    dS = ((-beta * S * I) / N)
    dI = ((beta * S * I) / N) - (gamma * I)
    dR = (gamma * I)
    dJ = ((beta * S * I) / N)
    return dS, dI, dR, dJ

# Initial conditions are S0, I0, R0
# Integrate the SIR equations over the time grid, t.
solve = odeint(deriv, (S0, I0, R0, J0), t, args=(N, beta, gamma))
S, I, R, J = solve.T

如您所见，beta 我留空了，注释掉了。如果我有其他所有的值，并且知道在流行高峰期，10%的人口被感染，是否可以从所有信息中找到beta？我试过的是：

sol= solve_ivp(lambda beta: deriv,
               [t], t_eval= t)
print(sol)

但是我发现这种语法不起作用。我的方法有什么问题？我如何估算 beta 值？

Answer 1

这里最简单的方法是通过 beta 参数化上面的代码并绘制结果，即作为 beta 函数的峰值感染，然后查看它越过阈值的位置。定义函数：

def peak_infections_pct(beta, n_days_total = 100):

    # Total population, N.
    N = 1000
    # Initial number of infected and recovered individuals, I0 and R0.
    I0, R0 = 10, 0
    # Everyone else, S0, is susceptible to infection initially.
    S0 = N - I0 - R0
    J0 = I0
    # Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
    gamma = 1/7
    # A grid of time points (in days)
    t = np.linspace(0, n_days_total, n_days_total+1)

    # The SIR model differential equations.
    def deriv(y, t, N, beta, gamma):
        S, I, R, J = y
        dS = ((-beta * S * I) / N)
        dI = ((beta * S * I) / N) - (gamma * I)
        dR = (gamma * I)
        dJ = ((beta * S * I) / N)
        return dS, dI, dR, dJ

    # Initial conditions are S0, I0, R0
    # Integrate the SIR equations over the time grid, t.
    solve = odeint(deriv, (S0, I0, R0, J0), t, args=(N, beta, gamma))
    S, I, R, J = solve.T

    return np.max(I)/N

计算并绘图：

betas = np.linspace(0,1,101,endpoint = True)
peak_inf = [peak_infections_pct(b) for b in betas]
plt.plot(betas, peak_inf)
plt.plot(betas, 0.1*np.ones(len(betas)))

获得

所以答案大约是 beta ~ 0.25 更准确地说，只需求解 beta：

from  scipy.optimize import root
root(lambda b: peak_infections_pct(b)-0.1, x0 = 0.5).x

输出：

array([0.23847079])

请注意，我将时间间隔留作函数的输入——您可能需要使用不同的长度，因为流行病可能会持续 100 天以上

为了仔细检查，让我们将感染绘制为我们的 beta=0.2384.. 的时间函数：

确实峰值在 100（是 10%）

在给定其他条件的情况下，在特定时间点估计 ODE 中的参数

Estimating a parameter in an ODE at a certain time point, given other conditions

python

scipy

ode

solver