如何将自定义ODE模型接入EasyABC进行ABC-SMC参数估计？

我正好做过类似的案例，把你的ODE模型接入EasyABC其实只需要对齐ABC-SMC的核心要求：一个能根据输入参数生成与观测数据维度匹配的模拟结果的函数，以及参数的先验分布定义。下面是完整的可运行示例，基于你的代码修改优化：

1. 加载必要的包 #

首先确保安装并加载deSolve和EasyABC：

# 安装并加载依赖包
if (!"deSolve" %in% rownames(installed.packages())) {
  install.packages("deSolve")
}
if (!"EasyABC" %in% rownames(installed.packages())) {
  install.packages("EasyABC")
}
library(deSolve)
library(EasyABC)

2. 修正ODE模型与模拟函数 #

你的原模拟函数用比例取时间点容易出现偏差，直接指定实验观测的时间点更可靠，同时要包含0h的初始值，保证和实验数据维度完全对齐：

# ODE核心方程（和你的原代码一致）
ode_model <- function(time, init, params) {
  with(as.list(c(init, params)), {
    dContamination <- Contamination * r * (1 - Contamination/C) - d * exp(-g * time) * Contamination
    return(list(dContamination))
  })
}

# 适配EasyABC的模拟函数：输入参数，返回对应观测时间点的细菌数量
simulate_data <- function(params) {
  # 初始污染量（对应0h的观测值）
  init <- c(Contamination = 59)
  # 实验观测的时间点：0h,1h,2h,4h,8h,24h
  obs_times <- c(0, 1, 2, 4, 8, 24)
  # 求解ODE
  sim <- ode(y = init, times = obs_times, func = ode_model, parms = params, 
             method = "radau", atol = 1e-4, rtol = 1e-4)
  # 返回模拟的细菌数量（仅保留数值列，去掉时间列）
  return(sim[, 2])
}

3. 准备实验观测数据 #

替换成你自己的真实实验数据即可：

# 示例实验数据：对应0h,1h,2h,4h,8h,24h的细菌数量
obs_data <- c(59, 32, 28, 22, 19, 26)

4. 定义参数的先验分布 #

EasyABC要求每个参数的先验是一个返回单个样本的函数，和你之前的试取值范围保持一致：

prior_list <- list(
  r = function() runif(1, 0.1, 4),    # 你的trial_r范围
  C = function() runif(1, 6, 15),     # 你的trial_C范围
  d = function() runif(1, 1, 4),      # 你的trial_d范围
  g = function() runif(1, 0.01, 1)    # 你的trial_g范围
)
# 注：你原代码中的`trial_l`参数未在ODE模型中使用，故暂不加入；若需使用，只需在列表中添加对应先验并更新ODE方程即可

5. 运行ABC-SMC算法 #

调用ABC_smc函数，核心参数说明：

• model：你的模拟函数（输入参数，输出模拟数据）
• prior：先验分布列表
• nb_simul：每个SMC迭代的模拟次数（根据计算资源调整）
• n_particles：最终得到的后验样本数量
• summary_stat_target：实验观测数据（作为匹配目标）

# 运行ABC-SMC
abc_result <- ABC_smc(
  model = simulate_data,
  prior = prior_list,
  nb_simul = 5000,  # 每个迭代的模拟次数，算力充足可适当调高
  n_particles = 200, # 后验样本数量，数量越多结果越稳定
  summary_stat_target = obs_data,
  tolerance = c(0.2, 0.1, 0.05), # SMC迭代的容忍度阈值，逐步降低逼近后验
  verbose = TRUE
)

6. 分析与可视化结果 #

运行完成后，可查看参数的后验分布：

# 查看后验样本的统计信息
summary(abc_result$param)

# 绘制参数后验分布直方图
par(mfrow = c(2, 2))
hist(abc_result$param[, "r"], main = "Posterior Distribution of r", xlab = "r")
hist(abc_result$param[, "C"], main = "Posterior Distribution of C", xlab = "C")
hist(abc_result$param[, "d"], main = "Posterior Distribution of d", xlab = "d")
hist(abc_result$param[, "g"], main = "Posterior Distribution of g", xlab = "g")

# 用ggplot2绘制更美观的分面图（可选）
library(ggplot2)
library(tidyr)
param_df <- as.data.frame(abc_result$param) %>%
  pivot_longer(cols = everything(), names_to = "Parameter", values_to = "Value")
ggplot(param_df, aes(x = Value)) +
  geom_histogram(bins = 30, fill = "steelblue", alpha = 0.7) +
  facet_wrap(~Parameter, scales = "free") +
  theme_minimal() +
  labs(title = "Posterior Distributions of ODE Parameters")

关键注意事项 #

1. 模拟函数的一致性：必须保证模拟函数返回的结果长度和观测数据完全一致，否则EasyABC无法计算距离。
2. 先验顺序匹配：先验列表的参数顺序要和模拟函数接收的参数顺序完全对应。
3. 容忍度调整：tolerance参数控制SMC迭代的严格程度，逐步降低的容忍度会让后验样本越来越接近真实分布，可根据结果调整迭代次数和阈值。

内容的提问来源于stack exchange，提问作者HCAI