Замените несколько циклов for с помощью функции

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

N<-5 # number of sites
sites<-LETTERS[seq(from=1,to=N)]
to.r<-rbind(sites)

p.move.r<-seq.int(0.05,0.95,by=0.1) # prob of moving to a new site
p.leave<-0.01*p.move.r # prob of leaving the system w/out returning
p.move.out<-0.01*p.move.r # prob of moving in/out
p.stay<-1-(p.move.r+p.leave+p.move.out) # prob of staying in the same site 

Для этого примера я включил только 50 симуляций, но на самом деле я хотел бы иметь по крайней мере 1000 симуляций,

set.seed(13973)

reps<-50 # number of replicates/simulations
steps<-100 # number of time steps (hours, days, weeks, etc)
random<-runif(10000,0,1) # generating numbers from a random distribution

# Construct empty df to fill with data

rep.movements<-matrix(NA,nrow=reps,ncol=steps)
colnames(rep.movements)<-c(1:steps);rownames(rep.movements)<-c(1:reps)

rep.use<-matrix(NA,nrow=reps,ncol=N)
colnames(rep.use)<-c(reefs);rownames(rep.use)<-c(1:reps)

# Outer loop to run each of the initial parameters

for(w in 1:length(p.stay)){
     p.move<-matrix((p.move.r[w]/(N-1)),N,N)
     diag(p.move)<-0

# Construction of distance matrix
move<-matrix(c(0),nrow=(N+2),ncol=(N+2),dimnames=list(c(sites,"NA","left"),c(sites,"NA","left")))
from<-array(0,c((N+2),(N+2)),dimnames=list(c(sites,"NA","left"),c(sites,"NA","left")))
to<-array(0,c((N+2),(N+2)),dimnames=list(c(sites,"NA","left"),c(sites,"NA","left")))

# Filling movement-Matrix construction

for(from in 1:N){
    for(to in 1:N){
      if(from==to){move[from,to]<-p.stay[w]} else {move[from,to]<-p.move[from,to]}
      move[,(N+1)]<-(1-(p.leave[w]+p.move.out[w]))/N
      move[,(N+2)]<-(1-(p.leave[w]+p.move.out[w]))/N
      move[(N+1),]<-p.move.out[w]
      move[(N+2),]<-p.leave[w] 
}

}

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

 cumsum.move<-cumsum(data.frame(move)) # Cumulative sum of probabilities

В этой кумулятивной матрице буквы «A», «B», «C», «D» и «E» представляют разные сайты, «NA» представляет вероятность ухода и возвращения на будущем временном шаге, а «покинуло» представляет вероятность покинуть систему и не вернуться обратно. Затем я использую список случайных чисел для сравнения с совокупной матрицей вероятностей и определения «судьбы» этого конкретного животного.

for(o в 1:повторениях){

result<-matrix(as.character(""),steps) # Vector for storing sites
x<-sample(random,steps,replace=TRUE) # sample array of random number 
time.step<-data.frame(x) # time steps used in the simulation (i)
colnames(time.step)<-c("time.step")
time.step$event<-""

j<-sample(1:N,1,replace=T) # first column to be selected 
k<-sample(1:N,1,replace=T) # selection of column for ind. that move in/out  

for(i in 1:steps){
  for (t in 1:(N+1)){
    if(time.step$time.step[i]<cumsum.move[t,j]){
    time.step$event[i]<-to.r[t]
    break
   }
 }

 ifelse(time.step$event[i]=="",break,NA) 
 result[i]<-time.step$event[i]
 j<-which(to.r==result[i]) 
 if(length(j)==0){j<-k} 
}

result<-time.step$event

# calculate frequency/use for each replicate

use<-table(result)
use.tab<-data.frame(use)
use.tab1<-use.tab[-which(use.tab==""),]
mergeuse<-merge(use.tab2,use.tab,all.x=TRUE)
mergeuse[is.na(mergeuse)]<-0

# insert data into empty matrix

rep.movements[o,]<-result
rep.use[o,]<-mergeuse$Freq

}

} 
  # for the outer loop I have some matrices to store the results for each parameter,
  # but for this example this is not important
rep.movements
rep.use

Теперь основная проблема заключается в том, что для запуска всех симуляций для каждого начального параметра (в этом примере 10 значений) требуется много времени. Мне нужно найти лучший/более эффективный способ запустить 1000 симуляций/20 сайтов по всем исходным параметрам. Я не слишком знаком с функциями или другими способами ускорить эту задачу. Любые идеи или рекомендации будут оценены.

Заранее большое спасибо,


person user1626688    schedule 18.11.2012    source источник


Ответы (1)


Давайте сначала обернем ваш код в функцию. Я также добавил команды set.seed, чтобы сделать результат воспроизводимым. Вам нужно удалить их перед запуском симуляции.

sim1 <- function(reps=50, steps=100 ) {

  N<-5 # number of sites
  sites<-LETTERS[seq(from=1,to=N)]
  to.r<-rbind(sites)

  p.move.r<-seq.int(0.05,0.90,by=0.05) # prob of moving to a new site
  p.leave<-0.01*p.move.r # prob of leaving the system w/out returning
  p.move.out<-0.01*p.move.r # prob of moving in/out
  p.stay<-1-(p.move.r+p.leave+p.move.out) # prob of staying in the same site 

  set.seed(42)
  random<-runif(10000,0,1) # generating numbers from a random distribution

  cumsum.move <- read.table(text="A         B         C         D         E    NA.   left
                            A    0.0820000 0.3407822 0.6392209 0.3516242 0.3925942 0.1964 0.1964
                            B    0.1254937 0.4227822 0.6940040 0.3883348 0.4196630 0.3928 0.3928
                            C    0.7959865 0.8730183 0.7760040 0.7930623 0.8765180 0.5892 0.5892
                            D    0.8265574 0.8980259 0.8095507 0.8750623 0.9000000 0.7856 0.7856
                            E    0.9820000 0.9820000 0.9820000 0.9820000 0.9820000 0.9820 0.9820
                            NA.   0.9910000 0.9910000 0.9910000 0.9910000 0.9910000 0.9910 0.9910
                            left 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000 1.0000",header=TRUE)

  cumsum.move <- as.matrix(cumsum.move)

  for(o in 1:reps){

    result<-matrix(as.character(""),steps) # Vector for storing sites
    set.seed(42)
    x<-sample(random,steps,replace=TRUE) # sample array of random number 
    time.step<-data.frame(x) # time steps used in the simulation (i)
    colnames(time.step)<-c("time.step")
    time.step$event<-""

    set.seed(41)
    j<-sample(1:N,1,replace=T) # first column to be selected 
    set.seed(40)
    k<-sample(1:N,1,replace=T) # selection of column for ind. that move in/out  

    for(i in 1:steps){
      for (t in 1:(N+1)){
        if(time.step$time.step[i]<cumsum.move[t,j]){
          time.step$event[i]<-to.r[t]
          break
        }
      }

      ifelse(time.step$event[i]=="",break,NA) 
      result[i]<-time.step$event[i]
      j<-which(to.r==result[i]) 
      if(length(j)==0){j<-k} 
    }

    result<-time.step$event
  }
  result
}

Обратите внимание, что result перезаписывается во время каждой итерации поверх o. Я не думаю, что вы хотите этого, поэтому я исправил это. Кроме того, вы используете data.frame внутри цикла. Как правило, вы должны избегать data.frames внутри циклов, как чумы. Хоть они и очень удобны, но по эффективности ужасны.

sim2 <- function(reps=50, steps=100) {

  N<-5 # number of sites
  sites<-LETTERS[seq(from=1,to=N)]
  to.r<-rbind(sites)

  p.move.r<-seq.int(0.05,0.90,by=0.05) # prob of moving to a new site
  p.leave<-0.01*p.move.r # prob of leaving the system w/out returning
  p.move.out<-0.01*p.move.r # prob of moving in/out
  p.stay<-1-(p.move.r+p.leave+p.move.out) # prob of staying in the same site 

  set.seed(42)
  random<-runif(10000,0,1) # generating numbers from a random distribution

  cumsum.move <- read.table(text="A         B         C         D         E    NA.   left
                            A    0.0820000 0.3407822 0.6392209 0.3516242 0.3925942 0.1964 0.1964
                            B    0.1254937 0.4227822 0.6940040 0.3883348 0.4196630 0.3928 0.3928
                            C    0.7959865 0.8730183 0.7760040 0.7930623 0.8765180 0.5892 0.5892
                            D    0.8265574 0.8980259 0.8095507 0.8750623 0.9000000 0.7856 0.7856
                            E    0.9820000 0.9820000 0.9820000 0.9820000 0.9820000 0.9820 0.9820
                            NA.   0.9910000 0.9910000 0.9910000 0.9910000 0.9910000 0.9910 0.9910
                            left 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000 1.0000",header=TRUE)

  cumsum.move <- as.matrix(cumsum.move)

  res <- list()
  for(o in 1:reps){

    result<-character(steps) # Vector for storing sites
    set.seed(42)
    time.step<-sample(random,steps,replace=TRUE) # sample array of random number 
    #time.step<-data.frame(x) # time steps used in the simulation (i)
    #colnames(time.step)<-c("time.step")
    #time.step$event<-""
    event <- character(steps)

    set.seed(41)
    j<-sample(1:N,1,replace=T) # first column to be selected 
    set.seed(40)
    k<-sample(1:N,1,replace=T) # selection of column for ind. that move in/out  

    for(i in 1:steps){
      for (t in 1:(N+1)){
        if(time.step[i]<cumsum.move[t,j]){
          event[i]<-to.r[t]
          break
        }
      }

      ifelse(event[i]=="",break,NA) 
      result[i]<-event[i]
      j<-which(to.r==result[i]) 
      if(length(j)==0){j<-k} 
    }

    res[[o]]<-event
  }
  do.call("rbind",res)
}

Обе функции дают одинаковый результат?

res1 <- sim1()
res2 <- sim2()
all.equal(res1,res2[1,])
[1] TRUE

Новая версия быстрее?

library(microbenchmark)
microbenchmark(sim1(),sim2())

Unit: milliseconds
    expr       min        lq    median        uq       max
1 sim1() 204.46339 206.58508 208.38035 212.93363 269.41693
2 sim2()  77.55247  78.39698  79.30539  81.73413  86.84398

Ну, множитель три - это уже неплохо. Я не вижу больших возможностей для дальнейшего улучшения циклов из-за этих break. Это оставляет только распараллеливание в качестве опции.

sim3 <- function(ncore=1,reps=50, steps=100) {
  require(foreach)
  require(doParallel)


  N<-5 # number of sites
  sites<-LETTERS[seq(from=1,to=N)]
  to.r<-rbind(sites)

  p.move.r<-seq.int(0.05,0.90,by=0.05) # prob of moving to a new site
  p.leave<-0.01*p.move.r # prob of leaving the system w/out returning
  p.move.out<-0.01*p.move.r # prob of moving in/out
  p.stay<-1-(p.move.r+p.leave+p.move.out) # prob of staying in the same site 

  set.seed(42)
  random<-runif(10000,0,1) # generating numbers from a random distribution

  cumsum.move <- read.table(text="A         B         C         D         E    NA.   left
                            A    0.0820000 0.3407822 0.6392209 0.3516242 0.3925942 0.1964 0.1964
                            B    0.1254937 0.4227822 0.6940040 0.3883348 0.4196630 0.3928 0.3928
                            C    0.7959865 0.8730183 0.7760040 0.7930623 0.8765180 0.5892 0.5892
                            D    0.8265574 0.8980259 0.8095507 0.8750623 0.9000000 0.7856 0.7856
                            E    0.9820000 0.9820000 0.9820000 0.9820000 0.9820000 0.9820 0.9820
                            NA.   0.9910000 0.9910000 0.9910000 0.9910000 0.9910000 0.9910 0.9910
                            left 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000 1.0000",header=TRUE)

  cumsum.move <- as.matrix(cumsum.move)

  #res <- list()
  #for(o in 1:reps){
  cl <- makeCluster(ncore)
  registerDoParallel(cl)
  res <- foreach(1:reps) %dopar% {

    result<-character(steps) # Vector for storing sites
    set.seed(42)
    time.step<-sample(random,steps,replace=TRUE) # sample array of random number 
    #time.step<-data.frame(x) # time steps used in the simulation (i)
    #colnames(time.step)<-c("time.step")
    #time.step$event<-""
    event <- character(steps)

    set.seed(41)
    j<-sample(1:N,1,replace=T) # first column to be selected 
    set.seed(40)
    k<-sample(1:N,1,replace=T) # selection of column for ind. that move in/out  

    for(i in 1:steps){
      for (t in 1:(N+1)){
        if(time.step[i]<cumsum.move[t,j]){
          event[i]<-to.r[t]
          break
        }
      }

      ifelse(event[i]=="",break,NA) 
      result[i]<-event[i]
      j<-which(to.r==result[i]) 
      if(length(j)==0){j<-k} 
    }

    #res[[o]]<-event
    event
  }
  stopCluster(cl)
  do.call("rbind",res)
}

Тот же результат?

res3 <- sim3()
all.equal(res1,c(res3[1,]))
[1] TRUE

Быстрее? (Давайте воспользуемся 4 ядрами на моем Mac. Вы можете попытаться получить доступ к серверу с еще несколькими ядрами.)

microbenchmark(sim1(),sim2(),sim3(4))
Unit: milliseconds
     expr        min         lq     median         uq       max
1  sim1()  202.28200  207.64932  210.32582  212.69869  255.2732
2  sim2()   75.39295   78.95882   80.01607   81.49027  125.0866
3 sim3(4) 1031.02755 1046.41610 1052.72710 1061.74057 1091.2175

Это выглядит ужасно. Однако тест несправедлив по отношению к параллельной функции. Функция вызывается 100 раз только с 50 повторениями. Это означает, что мы получаем все накладные расходы на распараллеливание, но почти не получаем от этого никакой пользы. Давайте сделаем это более справедливо:

microbenchmark(sim1(rep=10000),sim2(rep=10000),sim3(ncore=4,rep=10000),times=1)
Unit: seconds
                          expr      min       lq   median       uq      max
1            sim1(rep = 10000) 42.16821 42.16821 42.16821 42.16821 42.16821
2            sim2(rep = 10000) 16.13822 16.13822 16.13822 16.13822 16.13822
3 sim3(ncore = 4, rep = 10000) 38.18873 38.18873 38.18873 38.18873 38.18873

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

person Roland    schedule 18.11.2012
comment
Большое спасибо, Роланд, я действительно внес некоторые изменения в исходный код, чтобы включить более подробный/воспроизводимый пример кода, который я использую (который также имеет некоторые дополнительные циклы...). Сначала я проверю ваш сценарий. Спасибо за ваш вклад! - person user1626688; 18.11.2012
comment
Привет, Роланд, я не уверен, что этот код будет полезен при включении более подробных правок... так как мне нужно запускать симуляции для каждого параметра. - person user1626688; 18.11.2012
comment
Я уверен, что ваш новый код можно оптимизировать аналогичным образом. Суть в том, чтобы использовать эффективные структуры данных (векторы или списки) для хранения значений. Используйте функцию Rprof, чтобы узнать, какая операция занимает больше всего времени, и попытайтесь оптимизировать этот код. По возможности используйте векторизованные функции. И посмотрите на распараллеливание. - person Roland; 18.11.2012
comment
Привет, Роланд, какова цель функции do.call(rbind,res)? - person user1626688; 19.11.2012
comment
Это объединяет все векторы в списке в матрицу или data.frame. - person Roland; 19.11.2012