如何从以下模型函数中删除循环?

How to remove loop from the following model function?

我正在重写一些代码,目前正在创建一个小人口模型。我从一本书中重新创建了下面的当前模型函数,它是一个基于几个参数的简单人口模型。我将它们保留为默认值并返回了数据框。一切正常。但是,我想知道我是否可以以某种方式从函数中排除循环。

我知道 R 很棒,因为它可以进行矢量化计算,但我不确定在这种情况下是否可行。我想过用 lead/lag 之类的东西来做,但这行得通吗?也许不是因为事情需要按顺序计算?



# Nt numbers at start of time t
# Ct = removed at the end of time t
# Nt0 = numbers at time 0
# r = intrinsic rate of population growth
# K = carrying capacity


mod_fun = function (r = 0.5, K = 1000, N0 = 50, Ct = 0, Yrs = 10, p = 1) 
{
  # sets years to year value plus 1
  yr1 <- Yrs + 1
  # creates sequence of length years from year 1 to Yrs value +!
  years <- seq(1, yr1, 1)
  # uses years length to create a vector of length Yrs + 1
  pop <- numeric(yr1)
  # sets population at time 0
  pop[1] <- N0
  
  # creates a loop that calculates  model for each year after first year 
  for (i in 2:yr1) {
    # sets starting value of population for step to one calculated previous step
    # thus Nt is always the previous step pop size
    Nt <- pop[i - 1]
    
    pop[i] <- max((Nt + (r * Nt/p) * (1 - (Nt/K)^p) - 
                     Ct), 0)
  }
  
  # sets pop2 to original pop length
  pop2 <- pop[2:yr1]
  
  # binds together years (sequence from 1 to length Yrs), 
  # pop which is created in loop and is the population at the start of step t
  # pop2 which is the population at the end of step t
  out <- cbind(year = years, nt = pop, nt1 = c(pop2, NA))
  
  # sets row names to 
  rownames(out) <- years
  
  out <- out[-yr1, ]
  
  #returns data.frame
  return(out)
}

  
result = mod_fun()

这是输出的样子。基本上从第 1 行开始,给定起始人口 50,循环计算 nt1,然后将下一个 nt 行设置为 lag(nt1),然后事情以类似的方式继续。

result
#>    year       nt      nt1
#> 1     1  50.0000  73.7500
#> 2     2  73.7500 107.9055
#> 3     3 107.9055 156.0364
#> 4     4 156.0364 221.8809
#> 5     5 221.8809 308.2058
#> 6     6 308.2058 414.8133
#> 7     7 414.8133 536.1849
#> 8     8 536.1849 660.5303
#> 9     9 660.5303 772.6453
#> 10   10 772.6453 860.4776

Created on 2022-04-24 by the reprex package (v2.0.1)

mod_fun = function (r = 0.5, K = 1000, N0 = 50, Ct = 0, Yrs = 10, p = 1) 
{
  years <- seq_len(Yrs)
  pop <- Reduce(function(Nt, y)max((Nt + (r * Nt/p) * (1 - (Nt/K)^p) - Ct), 0),
         years, init = N0, accumulate = TRUE)
  
  data.frame(year = years, nt = head(pop,-1), nt1 = pop[-1])
 
}

   year       nt      nt1
1     1  50.0000  73.7500
2     2  73.7500 107.9055
3     3 107.9055 156.0364
4     4 156.0364 221.8809
5     5 221.8809 308.2058
6     6 308.2058 414.8133
7     7 414.8133 536.1849
8     8 536.1849 660.5303
9     9 660.5303 772.6453
10   10 772.6453 860.4776