On 2012-08-27 22:25, David Winsemius wrote:
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.
Thanks for the suggestion.
Because in real situation the parameters of `rnorm' function is
different in every loop, so it can't using 10,000 by 1,000 matrix
outside the loop.
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]
A mistake here, the above line should be changed to:
for (i in 1:99) qmatrix[i,] <- sort(pmatrix[i,])[c(0.05,0.5,0.95)*10000]
Thanks again.
Regards,
Jinsong
David Winsemius, MD
Alameda, CA, USA
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
