图形近似的迭代

Question

下面的代码运行一个基于微分方程组解生成的图形的模型，并显示了向点添加随机 5% 误差的模型。我将如何迭代这个过程（比如 100 次）并收集每个随机图生成的数据？看来我每次都必须创建新模型，这是不可行的。

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

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
  ])

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], # used to be np.linspace(0, 60, 1000)[48]
          args=(b, s, c, p, d)
      )
      data = solution.y[2, 0] / 0.01950269536785707  
      solutions.append(data)
  return np.array(solutions)

b = 0.00001
s = 4
c = 4
p = 2000000
p0 = (b, s, c, p)

np.random.seed(20)
D = np.logspace(-3, 3, 7)
z = model(D, b, s, c, p)
s = np.ones_like(z) * 0.05
n = random.gauss(0,z*0.05)
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.98458777e-05, 3.39383754e+00, 4.55115392e+00, 1.00007348e+06]),
#  array([[ 8.35308599e-10, -3.25230641e-03,  3.73469971e-03,
#          -5.22169803e-11],
#         [-3.25230641e-03,  1.28672442e+04, -1.47634805e+04,
#           2.06436797e-04],
#         [ 3.73469971e-03, -1.47634805e+04,  1.69398903e+04,
#          -2.36868204e-04],
#         [-5.22169803e-11,  2.06436797e-04, -2.36868204e-04,
#           3.31209815e-12]]))

Dlog = np.logspace(-3, 3, 200)
fig, axe = plt.subplots()
axe.scatter(D, zn, label="Data")
axe.semilogx(Dlog, model(Dlog, *popt), label="Fit")
axe.semilogx(Dlog, model(Dlog, *p0), "--", label="Model")
axe.legend()
axe.grid()
plt.show()

Answer 1

因此，微分方程组由函数“剂量”定义：

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
])

这里我们有三个变量： ‘目标’、‘传染性’、‘病毒’

要迭代此过程 100 次并收集每个随机图生成的数据，需要循环执行整个建模和拟合过程，并存储每次迭代的结果。那就是：

首先，您需要定义一个函数来生成噪声数据并拟合模型。其次，创建一个循环来多次执行此过程{!important}并收集并存储每次迭代的结果。

为此，创建多次迭代：

    ## Parameters
b = 0.00001
s = 4
c = 4
p = 2000000
p0 = (b, s, c, p)
D = np.logspace(-3, 3, 7)
sigma = 0.05
iterations = 100

# Storage of results
results = []

for i in range(iterations):
    zn, popt, pcov = run_iteration(D, b, s, c, p, p0, sigma, seed=i)
    results.append((zn, popt

这是可以提供帮助的循环迭代：

for i in range(iterations):
    zn, popt, pcov = run_iteration(D, b, s, c, p, p0, sigma, seed=i)
    results.append((zn, popt, pcov))

模型模拟： 求解微分方程组以获得理论数据 (z)。 将高斯噪声添加到数据中以模拟实验测量 (zn)。

曲线拟合： 使用 curve_fit 将模型拟合到噪声数据并获得参数估计值 (popt)。

重复和数据收集： 多次重复模拟和曲线拟合过程（100 次迭代）。 收集并存储调整后的参数和噪声数据以供以后分析。

每次迭代都会向数据添加一个小的随机变化（噪声），从而可以分析模型拟合如何随噪声数据变化。这很有用。

当然，如果代码中的所有内容都正确的话。这只是一个想法，不是 100% 求解器，只是想法。