基于计算衰减率的化学品降解曲线填充R函数问题排查

原始问题 #

尝试用计算好的动态衰减率（0-7天、7-13天各对应特定值）填充化学品浓度降解曲线，现有时间和初始浓度数据，但代码插值逻辑存在问题，无法正确生成每日浓度序列。

原始代码 #

interpolate_ts <- function(dates, concentrations) {
  data <- data.frame(date = dates, conc = concentrations)
  data <- data[order(data$date), ]
  data$decay <- c(NA, diff(log(data$conc)) / diff(data$date))
  seq_dates <- seq(min(data$date), max(data$date), by = 1)
  interp_conc <- rep(NA, nrow(data))

  # 为每个日期插值浓度，若为输入日期则使用真实值，否则计算得出，我怀疑问题出在这里
  for (i in seq_along(seq_dates)) {
    if (seq_dates[i] %in% data$date) {
      interp_conc[i] <- data$conc[data$date == seq_dates[i]]
    } else {
      prev_date <- lag(seq_dates, n = 1, default = NA)
      prev_conc <- lag(interp_conc, n = 1, default = NA)
      decay <- data$decay[data$date == prev_date]
      interp_conc[i] <- prev_conc + (decay * prev_conc)
    }
  }
return(interp_conc)
}

times = c(0, 7, 13, 14, 21, 56)
obs_conc = c(1.9659406, 1.3826208, 0.8696116, 2.5825205, 1.6871537, 0.4838216)
ts_interp = interpolate_ts(times, obs_conc)

核心问题分析 #

1. 向量长度不匹配：interp_conc初始化为和原始数据行数一致的向量，但seq_dates是每日序列（长度57），直接赋值会导致索引越界。
2. 滞后值获取错误：循环中用lag(seq_dates, n=1)返回的是整个序列的滞后向量，而非当前步的前一个日期，应该直接用seq_dates[i-1]。
3. 衰减率匹配逻辑错误：data$decay存储的是相邻原始数据点之间的区间衰减率，不是单个日期的衰减率，需要先为每个日期区间绑定对应的衰减率。
4. 降解公式错误：你用对数计算的衰减率对应指数降解模型，正确公式应为C(t) = C(t-1) * exp(decay_rate)，而非线性累加的prev_conc + decay*prev_conc。
5. 匹配结果不确定性：data$date == seq_dates[i]可能返回多个匹配值，需用match确保只取唯一对应值。

修正后的代码 #

interpolate_ts <- function(dates, concentrations) {
  # 整理原始数据并排序
  data <- data.frame(date = dates, conc = concentrations)
  data <- data[order(data$date), ]
  
  # 计算相邻区间的衰减率：data$decay[i]代表从date[i]到date[i+1]的日衰减率
  data$decay <- c(diff(log(data$conc)) / diff(data$date), NA)
  
  # 生成完整的每日日期序列
  seq_dates <- seq(min(data$date), max(data$date), by = 1)
  interp_conc <- rep(NA, length(seq_dates))
  
  # 先填充已知的观测浓度
  known_indices <- match(data$date, seq_dates)
  interp_conc[known_indices] <- data$conc
  
  # 遍历每日序列，填充插值浓度
  for (i in 2:length(seq_dates)) {
    if (is.na(interp_conc[i])) {
      # 找到当前日期所属的原始数据区间
      interval_idx <- max(which(data$date <= seq_dates[i]))
      # 获取该区间对应的衰减率
      decay_rate <- data$decay[interval_idx]
      # 用指数衰减公式计算当前浓度
      interp_conc[i] <- interp_conc[i-1] * exp(decay_rate)
    }
  }
  
  # 返回包含日期和插值浓度的数据框（更直观）
  return(data.frame(date = seq_dates, conc = interp_conc))
}

# 测试代码
times = c(0, 7, 13, 14, 21, 56)
obs_conc = c(1.9659406, 1.3826208, 0.8696116, 2.5825205, 1.6871537, 0.4838216)
ts_interp = interpolate_ts(times, obs_conc)

# 查看结果
head(ts_interp, 10)

修正说明 #

• 调整decay的计算逻辑：data$decay[i]明确对应date[i]到date[i+1]的日衰减率，最后一个值设为NA避免越界。
• 初始化interp_conc长度与seq_dates一致，避免索引错误。
• 用match快速定位已知浓度的位置，填充更高效。
• 循环中通过which(data$date <= seq_dates[i])找到当前日期所属的原始区间，匹配对应的衰减率。
• 使用指数衰减公式，和之前对数计算衰减率的逻辑保持一致。

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