Векторизация симуляции

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

#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))
}
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], можно ли векторизовать его с помощью функции применения?

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


r vectorization
person ECII    schedule 26.05.2012    source источник
comment
apply функции - это в основном синтаксический сахар вокруг for циклов, вы, вероятно, не выиграете большую скорость таким образом. Вы можете попробовать пакет compiler или Rcpp.   -  person baptiste    schedule 27.05.2012

Ответы (1)


arrow_upward
5
arrow_downward

Трудно «векторизировать» итерационные вычисления, но это симуляция, а симуляции, вероятно, будут выполняться много раз. Итак, напишите это, чтобы выполнить все симуляции одновременно, добавив аргумент 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)
}
person Martin Morgan    schedule 26.05.2012