Исправление функции, выполняемой для x: y вместо 1: y - PullRequest
Новая
       9

Исправление функции, выполняемой для x: y вместо 1: y

1 голос
/ 14 марта 2020

Я определил функцию для вычисления отношения между высотой (h) и диаметром (dbh) деревьев на основе уравнений, извлеченных из 2 публикаций. Моя цель состоит в том, чтобы использовать отношения, установленные в статье 1 (Сянтао), для прогнозирования значений переменных в уравнении в статье 2 (Маршо и Чаве). Я хотел бы проверить, в каком диапазоне диаметров [x:y] сгенерированная nls() кривая бумаги 2 подходит для бумаги 1. В настоящее время я получаю сообщение об ошибке (я верю в plot())

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

, если я использую что-либо кроме x = 1 для [x:y] т.е. 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 дней. Спасибо за любую помощь.

1 Ответ

1 голос
/ 15 марта 2020

Ваша проблема заключается в этом разделе: 

  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
  }

Что происходит на первый цикл л oop? Ну, tmp_h - это длина 20, но в первом цикле dbh - это 80, и вы присваиваете номер tmp_h[dbh], который равен tmp_h[80]. К тому времени, когда l oop закончится, в tmp_h будут храниться правильные значения, но они будут в индексах 80:100. Так что tmp_h будет иметь числа 80: 100, сохраненные в первых 21 индексах, затем группу 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
  }

и это будет работать.

Однако здесь вам вообще не нужен al oop, поскольку 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)))

, вы получаете это:

enter image description here

...