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

Question

Итак, в то время как lag и lead в dplyr великолепны, я хочу имитировать тайм-серверы чего-то типа роста населения. Моя старая школа код будет выглядеть примерно так:

tdf <- data.frame(time=1:5, pop=50) 
for(i in 2:5){ 
    tdf$pop[i] = 1.1*tdf$pop[i-1] 
} 

, который производит

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

Я чувствую, что должен быть dplyr или tidyverse способ сделать это (так же, как я люблю свою цикл) ,

Но что-то вроде

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

, который был бы моей первой догадкой только производит

time pop 
1 1 NA 
2 2 55 
3 3 55 
4 4 55 
5 5 55 

Я чувствую, что я что-то очевидное .... что это не хватаю?

Примечание. Это тривиальный пример. В моих реальных примерах используются несколько параметров, многие из которых меняются во времени (я моделирую прогнозы в разных сценариях GCM), поэтому tidyverse оказывается мощным инструментом в приведя мои симуляции вместе.

Answer 1

Reduce (или его purrr варианты, если вы хотите), что вы хотите для кумулятивных функций, которые уже не имеют 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 

или с purrr,

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 

Answer 2

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

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

Или,

base = 50 

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

Answer 3

Что о функциях карты, т.е.

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

Answer 4

Проблема здесь состоит в том, что 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)) 

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

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

Answer 5

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

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)))