图近似的迭代

下面的代码运行一个基于微分方程组解生成的图形的模型，并显示了在点上添加了随机 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()

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

# We will store the fitting parameters and generated data in these lists
fit_parameters = []
generated_data = []

# Number of iterations for the process
num_iterations = 100

for i in range(num_iterations):
  # Generate noisy data
  n = random.gauss(0, z * 0.05)
  zn = z + n

  # Fit the model to the noisy data
  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
  )

  # Store the fitting parameters and generated data
  fit_parameters.append(popt)
  generated_data.append(zn)

# Convert the lists to numpy arrays for easier analysis
fit_parameters = np.array(fit_parameters)
generated_data = np.array(generated_data)

# --- Visualization (optional) ---
Dlog = np.logspace(-3, 3, 200)

# Plot all the generated data and their fits
fig, axe = plt.subplots()
for i in range(num_iterations):
    axe.semilogx(Dlog, model(Dlog, *fit_parameters[i]), color='grey', alpha=0.2)
    axe.scatter(D, generated_data[i], color='blue', alpha=0.2, s=5) 

# Plot the original model and data for comparison
axe.semilogx(Dlog, model(Dlog, *p0), "--", label="Model")
axe.scatter(D, z, label="Original Data", color='red')

axe.legend()
axe.grid()
plt.show()

This code will:

1. Iterate the process of adding noise and fitting the model num_iterations times (set to 100 in the example).
2. Store the fitting parameters ( popt ) and the generated noisy data ( zn ) for each iteration in the fit_parameters and generated_data lists, respectively.
3. Convert these lists into NumPy arrays for easier analysis after the loop.

The optional visualization part will then plot all the generated data points and their corresponding fits, alongside the original model and data for comparison. You can adjust the plotting parameters like color , alpha , and s to change the appearance of the lines and markers.

This approach avoids creating a new model each time and provides you with the fitted parameters and generated data for further analysis.