ускорить симуляцию Монте-Карло с помощью вложенного цикла в R

Я хотел бы ускорить приведенное ниже моделирование методом Монте-Карло оценки DEA.

A<-nrow(banks)
effm<-matrix(nrow=A, ncol=2)
m<-20
B<-100

pb <- txtProgressBar(min = 0,
                     max = A, style=3)
for(a in 1:A) {
  x1<-x[-a,]
  y1<-y[-a,]
  theta=matrix(nrow=B,ncol=1) 

  for(i in 1:B){

    xrefm<-x1[sample(1:nrow(x1),m,replace=TRUE),]
    yrefm<-y1[sample(1:nrow(y1),m,replace=TRUE),]
    theta[i,]<-dea(matrix(x[a,],ncol=3),
                   matrix(y[a,],ncol=3),
                   RTS='vrs',ORIENTATION='graph',
                   xrefm,yrefm,FAST=TRUE)
  }

  effm[a,1]=mean(theta)
  effm[a,2]=apply(theta,2,sd)/sqrt(B)
  setTxtProgressBar(pb, a) 
}
close(pb)
effm

Как только A становится большим, симуляция останавливается. Я знаю из онлайн-исследований, что функция применения быстро ускоряет такой код, но я не знаю, как использовать ее в описанной выше процедуре.

Любая помощь/направление будет высоко оценена

Барри


r nested-loops montecarlo
person barryq    schedule 29.11.2012    source источник
comment
В сети много дезинформации. Функция apply может работать быстрее, чем цикл for, а может и нет; это зависит от того, что вы делаете. Вам нужно профилировать свой код на скорость, чтобы увидеть, какие части самые медленные (см. ?Rprof), тогда вы будете знать, что должно быть быстрее. Люди могут помочь профилировать ваш код, если вы предоставите воспроизводимый пример.   -  person Joshua Ulrich    schedule 29.11.2012
comment
@JoshuaUlrich то же самое! Кроме того, если вы можете опубликовать часть данных, которые вы используете, мы сможем фактически запустить ваш код, что значительно упростит помощь.   -  person Justin    schedule 29.11.2012
comment
Можете ли вы определить заморозку? Есть большая разница между процессом, который занимает много времени, и процессом, который сбрасывает системную память (или что-то еще) и зависает процесс и/или вся ОС.   -  person Carl Witthoft    schedule 29.11.2012
comment
Было бы полезно, если бы мы могли запускать этот код локально. Что такое banks?   -  person Roman Luštrik    schedule 30.11.2012

Ответы (1)


arrow_upward
1
arrow_downward

Следующее должно быть быстрее .... но если вы блокируете, когда A велико, это может быть проблемой памяти, а следующее требует больше памяти. Было бы полезно получить дополнительную информацию, например, что такое banks, что такое x, y, откуда вы получаете dea и какова цель.

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

A <- nrow(banks)
effm <- matrix(nrow = A, ncol = 2)
m <- 20
B <- 100
pb <- txtProgressBar(min = 0,
                     max = A, style=3)
for(a in 1:A) {
  x1 <- x[-a,]
  y1 <- y[-a,]
  theta <- numeric(B)
  xrefm <- x1[sample(1:nrow(x1), m * B, replace=TRUE),] # get all of your samples at once
  yrefm <- y1[sample(1:nrow(y1), m * B, replace=TRUE),]
  deaX <- matrix(x[a,], ncol=3)
  deaY <- matrix(y[a,], ncol=3)

  for(i in 1:B){
    theta[i] <- dea(deaX, deaY, RTS = 'vrs', ORIENTATION = 'graph',
                   xrefm[(1:m) + (i-1) * m,], yrefm[(1:m) + (i-1) * m,], FAST=TRUE)
  }

  effm[a,1] <- mean(theta)
  effm[a,2] <- sd(theta) / sqrt(B)
  setTxtProgressBar(pb, a) 
}
close(pb)
effm
person John    schedule 29.11.2012
comment
Было бы еще быстрее, если бы вы создали первые два аргумента для dea перед вторым циклом. - person Joshua Ulrich; 29.11.2012
comment
Если на то пошло, xynrow<-nrow(x1);mB<-m*B и переписывание как xrefm <- x1[sample(1:xynrow, mB, replace=TRUE),] удалит из вашей задачи огромное количество вычислений. - person Carl Witthoft; 29.11.2012
comment
@John Джон, я думаю, что dea из пакета бенчмаркинга, но можете ли вы объяснить в двух предложениях свою идею выполнить код? - person agstudy; 29.11.2012
comment
Спасибо за все полезные комментарии и извинения за некоторую неопределенность в моем сообщении. Мой код выполняет оценку гиперболической эффективности порядка M, где - person barryq; 29.11.2012
comment
Спасибо за все полезные комментарии/предложения по коду и извинения за некоторую неопределенность в моем посте. Банки - это кадр данных из 1394 наблюдений банка, x - матрица входных данных 1394 x 3, а y - матрица выходных данных 1394, а DEA - из бенчмаркинга. Я пытаюсь запустить гиперболическую оценку эффективности порядка M, где анализируемая единица (в данном случае каждый банк) не входит в эталонный набор в оценке DEA. - person barryq; 29.11.2012
comment
Что я имею в виду под зависанием, так это то, что код зависает в R, но влияет на вычислительную мощность ОС. - person barryq; 29.11.2012
comment
Я попытаюсь создать воспроизводимую копию процесса, и, надеюсь, это облегчит предоставление советов ... Должен ли я тогда создать новый вопрос или добавить в этот поток дополнительный комментарий? - person barryq; 29.11.2012
comment
@john Большое спасибо за изменение кода, но, к сожалению, это выдает ошибку Ошибка в dea(deaArg1, deaArg2, RTS = vrs, ORIENTATION = graph, xrefm[(1:m) +: количество входов должно быть одинаковым в X и XREF - person barryq; 29.11.2012
comment
трудно отлаживать без x... Вы можете удалять комментарии, которые не можете редактировать. Вы должны помещать много этого в свой вопрос и удалять комментарии. - person John; 29.11.2012
comment
@john yes теперь работает, но R все еще зависает примерно на 31% процесса - person barryq; 30.11.2012