r - 修复要为x :y instead of only 1:y运行的函数

Question

我定义了一个函数，该函数根据从2种出版物中提取的方程式来计算树木的高度(h)和直径(dbh)之间的关系。我的目标是使用论文1(Xiangtao)中建立的关系来预测论文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)))

你得到这个: