如何在R中计算多次清洁后牛舍粪便的累计残留量

问题描述 #

我需要评估牛舍每日清洁频率对地面粪便累积量的影响，但现有R代码仅能正确计算首次清洁后的粪便量，后续时段无法继承清洁状态。例如：

• 第7小时粪便量应为1676kg（无清洁累计值）
• 第8小时清洁后，正确值应为(1676 + 239.5) * 0.05 = 95.8kg
• 第9小时应为95.8 + 239.5 = 335.3kg
  但现有代码算出第8小时为2155.5kg，完全不符合预期。

问题根源 #

原代码核心错误在于未基于上一小时的实际粪便存量递推计算，而是直接使用无清洁状态下的累计值noclean来计算清洁后的量，导致清洁事件的影响无法被后续时段继承，每次清洁都错误地基于从头累计的总量计算。

解决方案 #

需要逐小时递推计算，每小时的粪便量必须依赖上一小时的结果：

1. 初始化第1小时的粪便量（无清洁时即为每小时产粪量）
2. 从第2小时开始，根据是否有清洁事件分别计算：
   • 无清洁：当前粪便量 = 上一小时粪便量 + 每小时产粪量
   • 有清洁：当前粪便量 = (上一小时粪便量 + 每小时产粪量) * 残留比例H

修正后代码 #

# 常量定义
totalfeceshour <- 239.5  # 250头牛每小时产粪量
H <- 0.05  # 刮除地面后残留粪便比例（原代码1-0.95，直接用0.05更清晰）

# 变量定义
hour <- 1:24
cleanbin <- c(0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1)  # 清洁事件发生在第8、16、24小时

# 初始化清洁后的粪便量向量
clean <- numeric(length(hour))
clean[1] <- totalfeceshour  # 第1小时无清洁，粪便量为每小时产粪量

# 逐小时递推计算
for(i in 2:length(hour)) {
  if(cleanbin[i] == 1) {
    # 有清洁：先加当前小时产粪，再乘以残留比例
    clean[i] <- (clean[i-1] + totalfeceshour) * H
  } else {
    # 无清洁：直接累加当前小时产粪
    clean[i] <- clean[i-1] + totalfeceshour
  }
}

# 合并结果为数据框（比矩阵更易读）
result <- data.frame(
  hour = hour,
  cleanbin = cleanbin,
  noclean = cumsum(rep(totalfeceshour, 24)),  # 无清洁累计值，用cumsum替代Reduce更简洁
  clean = clean
)

# 查看结果
print(result)

结果验证 #

运行代码后，关键时段的结果符合预期：

• 第8小时clean值为95.775kg（约95.8kg）
• 第9小时clean值为335.275kg（约335.3kg）

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