I would like to know if there is a more efficiency way to speed up below code. It uses a procedure where subsampling is required in the nested loop (which a previous answer helped to make more efficient). R has a tendency to hang when B=500, although computer OS isn't unduly affected. The goal is to run the below code with B=1000 and using larger m values(m=75, m=100, m=150)I have detailed the procedure in the code below and included a link to a reproducible data set.#Estimation of order m `leave one out' hyperbolic efficiency scores#The procedure sequentially works through each observation in `IOs' and #calculates a DEA order M efficiency score by leaving out the observation# under investigation in the DEA reference set# Step 1: Load the packages, create Inputs (x) and Outputs (y), choose # m(the order of the partial frontier or reference set to be used in DEA)# and B the number of monte carlo simulations of each m order DEA estimate# Step 2: For each observations a in x, x1 and y1#are create which 'leaves out' this observation. # Step 3: From these matrices subsamples (xref, yref) of size [m,] are # taken and used in DEA estimation.# Step 4: The DEA estimation uses the m subsample from step 3# as a reference set and evaluates the efficiency of the observation that# has been 'left out'#(thus the first two arguments in DEA are matrices of order [1,3] )# Step 5: Steps 3 and 4 are repeated B times to obtain B simulations of the # order m efficiency score and a mean and standard deviation are# calculated and placed in effm.# IOs data can be found here: https://dl.dropbox.com/u/1972975/IOs.txt# From IOs an Input matrix (x[1376,3]) and an Output matrix (y[1376,3])# are created.library(Benchmarking)x <- IOs[,1:3]y<-IOs[,4:6]A<-nrow(x)effm <- matrix(nrow = A, ncol = 2)m <- 50B <- 500pb <- txtProgressBar(min = 0,                     max = A, style=3)for(a in 1:A) {  x1 <- x[-a,]  y1 <- y[-a,]  theta <- numeric(B)  xynrow<-nrow(x1)  mB<-m*B   xrefm <- x1[sample(1:xynrow, mB, replace=TRUE),] # get all of your samples at once(https://stackoverflow.com/a/13629611/1176697)  yrefm <- y1[sample(1:xynrow, mB, replace=TRUE),]  deaX <- as.matrix(x[a,], ncol=3)  deaY <-as.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)
Speed up a monte carlo simulation with nested loop
r;simulation