함수를 사용하여 for 루프를 여러 개 바꾸기

Question

누군가이 예제를 빠르게 살펴보고이 문제에보다 효율적으로 접근 할 수 있기를 바랍니다. 특정 조건 세트를 따르는 동물이 사이트 간 이동하는 방식을 조사하기 위해 시뮬레이션을 실행하고 싶습니다. 나는 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"는 향후 시간 단계에서 이탈하고 다시 돌아올 확률을 나타내며 "왼쪽"은 시스템을 나가서 되돌아 가지 않을 확률을 나타냅니다. 그런 다음 난수 목록을 사용하여 누적 확률 행렬과 비교하고 특정 동물의 "운명"을 결정합니다. (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 개의 사이트를 실행하는보다 효율적이고 효율적인 방법을 찾아야합니다. 나는이 작업을 빠르게 할 수있는 기능이나 다른 방법에 익숙하지 않다. 어떤 아이디어 나 권장 사항도 환영 할 것입니다.

미리 감사드립니다.

Answer 1

함수에서 먼저 코드를 마무리합시다. 또한 결과를 재현 가능하게하기 위해 set.seed 명령을 추가했습니다. 시뮬레이션을 실행하기 전에 제거해야합니다. result O를 통해 각각의 반복 동안 덮어 것을

는

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 
} 

참고. 네가 그렇게하기를 원하지 않는다고 생각 했어. 또한 루프 내에 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 

글쎄요, 3 인자는 이미 아주 좋습니다. 루프를 더 향상시킬 수있는 가능성은별로 없습니다. 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 

더 빠름? (의 내 Mac에서 4 개 코어를 사용하자. 당신은 몇 가지 더 코어 서버에 접근을 시도 할 수도 있습니다.) 끔찍

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 

합니다. 그러나이 테스트는 병렬 기능에 불공평합니다. 이 함수는 50 회 반복만으로 100 회 호출됩니다. 즉, 병렬화의 오버 헤드가 발생하지만 거의 이점이 없습니다. 이제 더 공정한 만들어 보자 : 여전히

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 

더 나은,하지만 인상적 없습니다. 반복 횟수와 단계 수가 추가로 증가하면 병행 함수가 좋아 보이지만 필요한지 여부는 알 수 없습니다.