rreplicate

In R,Computing average from 1000 replications


I'm trying to do a test 1000 times, sum all of the answers and divide all of it by 1000.

I'm trying something like this but it is not working ,but it's only giving me the value for it doing it 1 time.

for(i in 1:1000) {
  x <- sample(1:500, 30, replace = F)

  mhat <- stimateM <- 2*mean(x) - 1
  y <- abs(mhat-500)
  value =+ y
}

Solution

  • Here's a more R-ish solution that doesn't use a loop:

    x <- replicate(1000, sample.int(500,30,replace=F))
    y <- abs(2*colMeans(x)-501)
    
    sumValue = sum(y)
    divValue = mean(y)
    

    replicate will give us the desired output as a matrix, and we can then use the vectorized function colMeans.

    This method will perform much better as your number of replicates increases:

    loopSoln <- function(N) {
      value = list()
      for(i in 1:N) {
        x=sample(1:500,30,replace=F)
        mhat=stimateM=2*mean(x)-1
        y=abs(mhat-500)
        value[i]=y
      }
      sumValue = sum(unlist(value))
      divValue = sumValue/N
    }
    
    replicateSoln <- function(N) {
      x <- replicate(1000, sample.int(500,30,replace=F))
      y <- abs(2*colMeans(x)-501)
      sumValue = sum(y)
      divValue = mean(y)
    }
    
    (ltimes <- sapply(c(1e2,1e3,1e4,1e5), function(N) system.time(loopSoln(N))[3]))
    ## elapsed elapsed elapsed elapsed 
    ##   0.002   0.014   0.158   2.009 
    
    (rtimes <- sapply(c(1e2,1e3,1e4,1e5), function(N) system.time(replicateSoln(N))[3]))
    ## elapsed elapsed elapsed elapsed 
    ##   0.007   0.011   0.007   0.010 
    
    plot(ltimes~I(2:5), type = 'l', col = 2, xlab = "log_10 number of replicates", ylab = "elapsed time (sec)")
    lines(rtimes~I(2:5), col = 4)
    legend("topleft", lty = 1, col = c(2,4), 
      legend = c("Loop", "replicate + rowMeans"))
    

    Benchmarks