修复一个函数以运行 x:y 而不是仅 1:y

Question

我已经定义了一个函数来根据从 2 个出版物中提取的方程计算树木的高度 (h) 和直径 (dbh) 之间的关系。我的目标是使用论文 1（向涛）中建立的关系来预测论文 2（Marechaux 和 Chave）中方程中变量的值。我想测试一下[x:y]生成的纸2的nls()曲线适合纸1的直径范围。目前，我不断收到错误消息（我相信plot()）

Error in xy.coords(x, y, xlabel, ylabel, log) : 'x' and 'y' lengths differ

如果我对[x:y] 使用除 x=1 以外的任何东西，即dbh.min:dbh.max

我的功能如下：

# Plant.Functional.Type constants...
Dsb1 <- 2.09
Dsb2 <- 0.54
Db1 <- 0.93
Db2 <- 0.84
BDb1 <- 2.66
BDb2 <- 0.48
Eb1 <- 1.41
Eb2 <- 0.65
# # # # # # # # # # # # # # # # # # # # # # # # # # #
Generate.curve <- function(b1, b2, dbh.min, dbh.max){
# calculate Xiangtao's allometry...
  tmp_h <- c(dbh.min:dbh.max)
  for (dbh in dbh.min:dbh.max)
  {
    h = b1*dbh^(b2)
    tmp_h[dbh] = h
  }
# plot to check curve
  plot(dbh.min:dbh.max, tmp_h)

# define secondary function for Marechaux and Chave allometry
  h_fxn <- function(hlim,dbh,ah){
    h = hlim * (dbh / (dbh + ah))
    return(h)
  }

# use nonlinear least squares model to solve for ah and hlim
  # set model inputs
  start.ah <- 1 
  start.hlim <- 5
  tmp_v <- cbind(dbh.min:dbh.max,tmp_h)

tmp.fit <- nls(tmp_h ~ h_fxn(hlim,dbh.min:dbh.max,ah), start = list(hlim = start.hlim, 
                ah = start.ah), algorithm = "port", upper = list(hlim = 75, ah = 99))  
# seems to be no way of extracting ah and hlim from tmp.fit via subset
# extract manually and then check fit with
  # lines(dbh.min:dbh.max, hlim * (dbh.min:dbh.max/(dbh.min:dbh.max + ah)))
  # for equation h = hlim * (dbh / (dbh + ah)) from Marechaux and Chave
return(tmp.fit)
}
# # # # # # # # # # # # # # # # # # # # # # # # # # #

这很适合

Generate.curve(Dsb1,Dsb2,1,100)
lines(1:100, 36.75 * (1:100/(1:100 + 52.51)))

但我也希望能够检查曲线拟合范围，例如 [80:100]。 我一直试图弄清楚为什么Generate.curve(Dsb1,Dsb2,80,100) 现在返回错误大约 3 天。感谢您的帮助。

Answer 1

你的问题出在这部分：

  tmp_h <- c(dbh.min:dbh.max)
  for (dbh in dbh.min:dbh.max)
  {
    h = b1*dbh^(b2)
    tmp_h[dbh] = h
  }

想想当您将dbh.min 设置为80 并将dbh.max 设置为100 时会发生什么：

  tmp_h <- 80:100
  for (dbh in 80:100)
  {
    h = b1*dbh^(b2)
    tmp_h[dbh] = h
  }

循环的第一个循环会发生什么？好吧，tmp_h 的长度为 20，但在第一个循环中，dbh 为 80，而您正在为tmp_h[dbh] 分配一个数字，即tmp_h[80]。到循环结束时，tmp_h 将存储正确的值，但它们将在索引80:100 中。所以tmp_h 将在前 21 个索引中存储数字 80:100，然后是一堆 NA，然后是最后 21 个索引中的正确数字。

所以改成：

  tmp_h <- c(dbh.min:dbh.max)
  for (dbh in dbh.min:dbh.max)
  {
    h = b1*dbh^(b2)
    tmp_h[dbh - dbh.min + 1] = h
  }

它会起作用的。

但是，这里实际上根本不需要循环，因为 R 使用矢量化操作，所以整个部分可以替换为：

tmp_h <- b1 * (dbh.min:dbh.max)^(b2)

然后当你这样做时

Generate.curve(Dsb1,Dsb2,80,100)
lines(80:100, 36.75 * (80:100/(80:100 + 52.51)))

你明白了：