[R] for loop optimization help

[David Winsemius]

On Aug 27, 2012, at 1:53 AM, Jinsong Zhao wrote:

> On 2012-08-27 9:35, David Winsemius wrote:
>>
>> On Aug 26, 2012, at 5:06 PM, Jinsong Zhao wrote:
>>
>>> Hi there,
>>>
>>> In my code, there is a for loop like the following:
>>>
>>>  pmatrix <- matrix(NA, nrow = 99, ncol = 10000)
>>>  qmatrix <- matrix(NA, nrow = 99, ncol = 3)
>>>  paf <- seq(0.01, 0.99, 0.01)
>>>  for (i in 1:10000) {
>>>      p.r.1 <- rnorm(1000, 1, 0.5)
>>>      p.r.2 <- rnorm(1000, 2, 1.5)
>>>      p.r.3 <- rnorm(1000, 3, 1)
>>>      pmatrix[,i] <- quantile(c(p.r.1, p.r.2, p.r.3), paf)
>>>  }
>>>  for (i in 1:99) {
>>>     qmatrix[i,] <- quantile(pmatrix[i,], c(0.05, 0.5, 0.95))
>>>  }
>>>
>>> Because of the number of loop is very large, e.g., 10000 here, the
>>> code is very slow.
>>
>> I would think that picking the seq(0.01, 0.99, 0.01) items in the  
>> first
>> case and the 500th, 5000th and the 9500th in the second case,  rather
>> than asking for what `quantile` would calculate,  would surely be  
>> more
>> "statistical", in the sense of choose order statistics anyway. Likely
>> much faster.
>>

Also possible that drawing the normals as 10,000 by 1,000 matrices  
(all at once, outside the loop) would save time (at the cost of added  
space requirements.

>
> Yes, you are right. Following your suggestions, the execution time  
> of `sort' is much shorter than `quantile' in the following code:
>
>  pmatrix <- matrix(NA, nrow = 99, ncol = 10000)
>  qmatrix <- matrix(NA, nrow = 99, ncol = 3)
>  paf <- seq(0.01, 0.99, 0.01)
>  for (i in 1:10000) {
>      p.r.1 <- rnorm(1000, 1, 0.5)
>      p.r.2 <- rnorm(1000, 2, 1.5)
>      p.r.3 <- rnorm(1000, 3, 1)
>      pmatrix[,i] <- sort(c(p.r.1, p.r.2, p.r.3))[paf*3000]
>  }
>  qmatrix <- pmatrix[,c(0.05, 0.5, 0.95)*10000]
>
> Thanks again.
>
> Regards,
> Jinsong

David Winsemius, MD
Alameda, CA, USA