如何修复用于估算模型年龄的R嵌套循环脚本

问题描述 #

我是R（及编程领域）的新手，所使用的数据来自铅同位素分析。我希望编写一个嵌套循环，用于对比某方程的输出与基于样本数据计算出的整数顶点。目前遇到两个问题：

1. 1:length(num_intervals)引用的是整数，无法生成正确的区间序列；
2. 循环的输出错误覆盖了ModelAge列的值。

修正后的完整代码 #

# 读取测试数据
ingots <- read.table(text = "
47      18.54400       15.590       38.427
48      18.55900       15.580       38.403
49      18.48100       15.600       38.407
50      18.37400       15.614       38.293
51      19.23800       15.673       37.272
52      18.14200       15.649       38.137
53      18.52500       15.590       38.430
54      18.55000       15.597       38.468
55      18.59200       15.662       37.972
56      18.59500       15.662       37.972
57      18.50000       15.594       38.402
", header = FALSE, col.names = c("SampleID", "X206Pb.204Pb", "X207Pb.204Pb", "X208Pb.204Pb"))

# Model Age Estimations
rm(list = ls())
# library(dplyr)

# 常数假设
decay_constant_u238 <- 1.55125E-10
decay_constant_u235 <- 9.8485E-10
decay_constant_th232 <- 4.9475E-11

# 时间假设
age_of_earth_t0 <- 4.57E9
age_of_earth_t1 <- 3.7E9

# 同位素假设 - Stage 1
s1_206PB.204PB <- 9.307
s1_207PB.204PB <- 10.294
s1_208PB.204PB <- 29.476

# 同位素假设 - Stage 2 (Stacey & Kramer 1975)
s2_206PB.204PB <- 11.152
s2_207PB.204PB <- 12.998
s2_208PB.204PB <- 31.230

# 同位素比例假设 - Stage 1 (Stacey & Kramer 1975)
s1_u <- 7.19
s1_w <- 32.21
s1_k <- 4.479833102

# 同位素比例假设 - Stage 2 (Stacey & Kramer 1975)
s2_u <- 9.74
s2_w <- 37.19
s2_k <- 3.818275154

# µ or 238U/204Pb Estimates
uranium_lead_ratio_estimates <- matrix(6:12, nrow=7, ncol=1)

# 计算斜率（补充原代码中缺失的变量定义）
values_207.204 <- ingots$X207Pb.204Pb
values_206.204 <- ingots$X206Pb.204Pb
slopes <- (values_207.204 - s2_207PB.204PB)/(values_206.204 - s2_206PB.204PB)*137.88

# 初始化结果数据框
model_age <- data.frame(
  A = ingots$X206Pb.204Pb, 
  B = ingots$X207Pb.204Pb, 
  C = ingots$X208Pb.204Pb, 
  Slopes = slopes, 
  ModelAge = numeric(length(ingots$X206Pb.204Pb)) # 初始化为数值型向量
)

# 函数定义
ModelAgeFunction <- function(Ages){
  (exp(decay_constant_u235*age_of_earth_t1)-exp(decay_constant_u235*Ages))/
    (exp(decay_constant_u238*age_of_earth_t1)-exp(decay_constant_u238*Ages))
}

# 迭代参数设置
lower_bound <- 0
upper_bound <- 700000000
num_intervals <- 70
interval_width <- (upper_bound - lower_bound) / num_intervals

# 嵌套循环修正
for (i in 1:length(slopes)) {
  # 遍历每个区间，修正为1:num_intervals
  for (j in 1:num_intervals) {
    # 基于内层循环的j计算区间左边界，而非外层的i
    left_bound <- lower_bound + interval_width*(j-1)
    mid_point <- left_bound + interval_width / 2
    if (slopes[i] < ModelAgeFunction(mid_point)) {
      # 用外层循环的i索引对应样本的ModelAge，而非j
      model_age$ModelAge[i] <- mid_point
      break
    }
  }
}

# 查看结果
print(model_age)

问题修复说明 #

1. 解决1:length(num_intervals)的序列问题 #

原代码中num_intervals是单个整数（70），length(num_intervals)返回1，导致内层循环仅执行一次。将内层循环的1:length(num_intervals)改为1:num_intervals，即可生成1到70的完整区间序列，实现对所有区间的遍历。

2. 解决输出覆盖问题 #

• 索引错误修正：原代码中用内层循环的j索引model_age$ModelAge[j]，会导致所有样本的结果都覆盖到同一个位置。改为用外层循环的i索引model_age$ModelAge[i]，确保每个样本的结果对应到数据框的正确行。
• 区间边界计算修正：原代码中left_bound基于外层循环的i计算，导致区间遍历错误。改为基于内层循环的j计算，才能正确遍历每个区间的左边界。
• 初始化修正：将ModelAge列初始化为数值型向量，避免原代码中用整数序列初始化可能导致的类型问题。

3. 补充原代码缺失部分 #

添加了测试数据的读取代码，以及values_207.204、values_206.204的变量定义，解决原代码运行时的未定义变量报错问题。

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