Моделирование временного ряда в dplyr вместо использования цикла for

Reduce (или его муррр-варианты, если хотите) — это то, что вам нужно для кумулятивных функций, для которых еще не написана версия cum*:

data.frame(time = 1:5, pop = 50) %>%
    mutate(pop = Reduce(function(x, y){x * 1.1}, pop, accumulate = TRUE))

##   time    pop
## 1    1 50.000
## 2    2 55.000
## 3    3 60.500
## 4    4 66.550
## 5    5 73.205

или с мурлыканьем,

data.frame(time = 1:5, pop = 50) %>%
    mutate(pop = accumulate(pop, ~.x * 1.1))

##   time    pop
## 1    1 50.000
## 2    2 55.000
## 3    3 60.500
## 4    4 66.550
## 5    5 73.205

Если начальное значение pop равно, скажем, 50, то pop = 50 * 1.1^(0:4) даст вам следующие четыре значения. С вашим кодом вы можете сделать:

data.frame(time=1:5, pop=50) %>%
  mutate(pop = pop * 1.1^(1:n() - 1))

Or,

base = 50

data.frame(time=1:5) %>%
  mutate(pop = base * 1.1^(1:n()-1))

Функция накопления Purrr может обрабатывать изменяющиеся во времени индексы, если вы передаете их в свою функцию моделирования в виде списка со всеми параметрами в нем. Тем не менее, нужно немного повозиться, чтобы заставить это работать правильно. Хитрость здесь в том, что функция calculate() может работать как со списком, так и с векторными столбцами. Вы можете использовать функцию nest() tidyr, чтобы сгруппировать столбцы в вектор списка, содержащий текущее состояние населения и параметры, а затем использовать функцию calculate() для результирующего столбца списка. Это немного сложно объяснить, поэтому я включил демонстрацию, моделирующую логистический рост либо с постоянной скоростью роста, либо с изменяющейся во времени стохастической скоростью роста. Я также включил пример того, как использовать это для имитации нескольких повторений для данной модели, используя dpyrr+purrr+tyr.

library(dplyr)
library(purrr)
library(ggplot2)
library(tidyr)

# Declare the population growth function. Note: the first two arguments
# have to be .x (the prior vector of populations and parameters) and .y,
# the current parameter value and population vector. 
# This example function is a Ricker population growth model. 
logistic_growth = function(.x, .y, growth, comp) {
  pop = .x$pop[1]
  growth = .y$growth[1]
  comp  = .y$comp[1]
  # Note: this uses the state from .x, and the parameter values from .y.
  # The first observation will use the first entry in the vector for .x and .y
  new_pop = pop*exp(growth - pop*comp)
  .y$pop[1] = new_pop
  return(.y)
}

# Starting parameters the number of time steps to simulate, initial population size,
# and ecological parameters (growth rate and intraspecific competition rate)
n_steps  = 100
pop_init = 1
growth = 0.5
comp = 0.05

#First test: fixed growth rates
test1 = data_frame(time = 1:n_steps,pop = pop_init, 
                   growth=growth,comp =comp)


# here, the combination of nest() and group_by() split the data into individual 
# time points and then groups all parameters into a new vector called state.
# ungroup() removes the grouping structure, then accumulate runs the function
#on the vector of states. Finally unnest transforms it all back to a
#data frame
out1 = test1 %>%
  group_by(time)%>%
  nest(pop, growth, comp,.key = state)%>%
  ungroup()%>%
  mutate(
    state = accumulate(state,logistic_growth))%>%
  unnest()

# This is the same example, except I drew the growth rates from a normal distribution
# with a mean equal to the mean growth rate and a std. dev. of 0.1
test2 = data_frame(time = 1:n_steps,pop = pop_init, 
                   growth=rnorm(n_steps, growth,0.1),comp=comp)

out2 = test2 %>%
  group_by(time)%>%
  nest(pop, growth, comp,.key = state)%>%
  ungroup()%>%
  mutate(
    state = accumulate(state,logistic_growth))%>%
  unnest()

# This demostrates how to use this approach to simulate replicates using dplyr
# Note the crossing function creates all combinations of its input values
test3 = crossing(rep = 1:10, time = 1:n_steps,pop = pop_init, comp=comp) %>%
  mutate(growth=rnorm(n_steps*10, growth,0.1))

out3 = test3 %>%
  group_by(rep)%>%
  group_by(rep,time)%>%
  nest(pop, growth, comp,.key = state)%>%
  group_by(rep)%>%
  mutate(
    state = accumulate(state,logistic_growth))%>%
  unnest()

print(qplot(time, pop, data=out1)+
        geom_line() +
        geom_point(data= out2, col="red")+
        geom_line(data=out2, col="red")+
        geom_point(data=out3, col="red", alpha=0.1)+
        geom_line(data=out3, col="red", alpha=0.1,aes(group=rep)))

Проблема здесь в том, что dplyr выполняет это как набор векторных операций, а не оценивает термин по одному. Здесь 1.1*lag(pop) интерпретируется как «вычислить запаздывающие значения для всего поп-музыки, а затем умножить их все на 1,1». Поскольку вы set pop=50 запаздывали, значения для всех шагов были 50.

dplyr имеет несколько вспомогательных функций для последовательной оценки; стандартные функции cumsum, cumprod и т. д. работают, а несколько новых (см. ?cummean) работают в dplyr. В вашем примере вы можете смоделировать модель с помощью:

tdf <- data.frame(time=1:5, pop=50, growth_rate = c(1, rep(1.1,times=4)) %>%
    mutate(pop = pop*cumprod(growth_rate))


time    pop     growth_rate
1       50.000  1.0
2       55.000  1.1
3       60.500  1.1
4       66.550  1.1
5       73.205  1.1

Обратите внимание, что здесь я добавил скорость роста в качестве столбца и установил для первой скорости роста значение 1. Вы также можете указать его следующим образом:

tdf <- data.frame(time=1:5, pop=50, growth_rate = 1.1) %>%
    mutate(pop = pop*cumprod(lead(growth_rate,default=1))

Это делает явным, что столбец скорости роста относится к скорости роста на текущем временном шаге по сравнению с предыдущим.

Существуют ограничения на количество различных симуляций, которые вы можете выполнить таким образом, но должно быть возможно построить множество экологических моделей с дискретным временем, используя некоторую комбинацию кумулятивных функций и параметров, указанных в столбцах.

Как насчет функций карты, т.е.

tdf <- data_frame(time=1:5)
tdf %>% mutate(pop = map_dbl(.x = tdf$time, .f = (function(x) 50*1.1^x)))