Векторизация моделирования

Question

Пытаясь обернуть вокруг вектора век моего разума, пытаясь сделать некоторые симуляции быстрее, я нашел это очень основное моделирование эпидемии. Код из книги http://www.amazon.com/Introduction-Scientific-Programming-Simulation-Using/dp/1420068725/ref=sr_1_1?ie=UTF8&qid=1338069156&sr=8-1

#program spuRs/resources/scripts/SIRsim.r 

SIRsim <- function(a, b, N, T) { 
    # Simulate an SIR epidemic 
    # a is infection rate, b is removal rate 
    # N initial susceptibles, 1 initial infected, simulation length T 
    # returns a matrix size (T+1)*3 with columns S, I, R respectively 
    S <- rep(0, T+1) 
    I <- rep(0, T+1) 
    R <- rep(0, T+1) 
    S[1] <- N 
    I[1] <- 1 
    R[1] <- 0 
    for (i in 1:T) { 
    S[i+1] <- rbinom(1, S[i], (1 - a)^I[i]) 
    R[i+1] <- R[i] + rbinom(1, I[i], b) 
    I[i+1] <- N + 1 - R[i+1] - S[i+1] 
    } 
    return(matrix(c(S, I, R), ncol = 3)) 
} 

Ядро моделирования является for цикла. Мой вопрос, так как код генерирует значения S[i+1] и R[i+1] из значений S[i] и R[i], можно ли его векторизовать с помощью функции приложения?

Большое спасибо

Answer 1

Это трудно «векторизация» итеративные вычисления, но это моделирование и моделирование, вероятно, будет работать во много раз. Поэтому напишите это, чтобы сделать все симуляции одновременно, добавив аргумент M (количество выполняемых имитаций), выделяя матрицу M x (T + 1), а затем заполняя последовательные столбцы (времена) каждого моделирования. Изменения выглядят удивительно прямолинейными (поэтому я, вероятно, допустил ошибку, я особенно обеспокоен использованием векторов во втором и третьем аргументах для rbinom, хотя это согласуется с документацией).

SIRsim <- function(a, b, N, T, M) { 
    ## Simulate an SIR epidemic 
    ## a is infection rate, b is removal rate 
    ## N initial susceptibles, 1 initial infected, simulation length T 
    ## M is the number of simulations to run 
    ## returns a list of S, I, R matricies, each M simulation 
    ## across T + 1 time points 
    S <- I <- R <- matrix(0, M, T + 1) 
    S[,1] <- N 
    I[,1] <- 1 
    for (i in seq_along(T)) { 
     S[,i+1] <- rbinom(M, S[,i], (1 - a)^I[,i]) 
     R[,i+1] <- R[,i] + rbinom(M, I[,i], b) 
     I[,i+1] <- N + 1 - R[,i+1] - S[,i+1] 
    } 
    list(S=S, I=I, R=R) 
}