计算参数的ARE

我有以下代码：

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

我想计算 ODE (b,s,c,p) 中参数的 ARE（平均残差）。我怎样才能有效地做到这一点？作为参考，这是通过以下计算得出的，其中 $\theta$ 是新的猜测，$\theta_i$ 是旧的猜测。

我想从图中读出参数：

当然，我可以帮助计算 ODE 参数 (b, s, c, p) 的平均相对误差 (ARE)。

理解 ARE

平均相对误差 (ARE) 衡量估计参数值与真实参数值之间的平均偏差百分比。其计算公式如下：

ARE = (1/n) * Σ [(|θ - θ̂|)/ |θ|] * 100%

其中：

• θ 表示真实的参数值
• θ̂ 表示估计的参数值
• n 表示参数的个数

在代码中计算 ARE

需要修改代码以包含 ARE 计算。以下是如何操作：

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

# ... （的现有代码） ...

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
)

# 计算 ARE
theta = np.array([b, s, c, p])  # 真实的参数值
theta_hat = popt                # 估计的参数值

are = (1/len(theta)) * np.sum(np.abs((theta - theta_hat) / theta)) * 100

print(f"平均相对误差 (ARE): {are:.2f}%")

# ... （代码的其余部分） ...

此代码段将计算并打印估计参数的 ARE。

从图中读取参数

提供的图像无法直接读取参数值。为了从图中获得参数的近似值，需要：

1. 确定曲线的函数形式: 观察数据点并尝试确定最适合数据的函数类型（例如，线性、指数、对数）。

2. 使用曲线拟合技术: 使用 Python 中的曲线拟合工具（例如， scipy.optimize.curve_fit 或 numpy.polyfit ）将选定的函数拟合到数据。这将为提供拟合参数的估计值。

重要说明: 从图中读取参数值本质上是不精确的。为了获得更准确的参数估计，需要使用原始数据和适当的拟合技术。