There are several things going on here.
First, you calculate the scales as:
scale.list <- runif(1000, max = 100000, min = 100000)
but since min = max
, all the values are identical.
Second, you do not specify lb.arg
or ub.arg
, so I set them to 20 and 50 arbitrarily.
Third, profiling this code using Rprof
shows that > 90% of the time is spent in the qtrunc(...)
function, which is called by rtrunc(...)
. This is because you are generating 100,000 sample points at each iteration, and qtrunc(...)
has to sort these. The total run time scales as O(n) where n is the number of sample points. On my system, using n=1000 takes about 7 sec, so using n=100,000 would take 700 sec or about 12 minutes.
My advice would be to try smaller n and see if that really makes a difference. From the central limit theorem we know that the distribution of the mean is asymptotically Normal for large n, regardless of the underlying distribution. I doubt that increasing n from 1000 to 100,000 changes that significantly.
Finally, the idiomatic way to do this in R is [using n=1000]:
f <- function(shape,scale,lb,ub){
X <- rtrunc(1000, spec="gamma", a=lb, b=ub,shape=shape,scale=scale)
return(c(mean=mean(X),sd=sd(X)))
}
# set.seed(1) # use this for a reproducible sample
result <- mapply(f,shape.list,scale.list, lb=20,ub=50)
result.apply <- data.frame(t(result))
which produces a data frame with two columns: the mean and the sd for each shape/scale. By setting the seed to a fixed value just prior to running mapply(...)
, and doing the same just prior to running your for loop, you can show that these both produce identical results.