向量化模擬

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]和R[i]值中的S[i+1]和R[i+1]值，是否可以使用apply函數對其進行矢量化？

非常感謝

Answer 1

很難「矢量化」迭代計算，但是這是一個模擬和仿真有可能被多次運行。因此，通過添加一個參數M（要執行的模擬次數），分配一個M x（T + 1）矩陣，然後填充每個模擬的連續列（時間），來編寫此代碼以同時執行所有模擬。這些變化看起來非常直截了當（所以我可能犯了一個錯誤;我特別關注在第二個和第三個參數中使用向量，但這與文檔一致）。

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