##
## A simple simulation to estimate mean and standard errors of the
## sample mean and sample median for a sample from a gamma
## distribution.
##

gsim <- function(R, n, a) {
    x <- matrix(rgamma(n * R, a), ncol = n)
    mn <- apply(x, 1, mean)
    md <- apply(x, 1, median)
    c(mn = mean(mn), md = mean(md), sd.mn = sd(mn), sd.md = sd(md))
}

## > system.time(val <- gsim(100000, 10, 2))
##    user  system elapsed 
##   5.320   0.004   5.336 

## a parallelized version that runs gsim on the nodes in the cluster
## and merges the results.
mergeMeans <- function(m, RR) sum(m * (RR / sum(RR)))
mergeSDs <- function(s, RR) sqrt(sum(s^2 * ((RR - 1) / (sum(RR) - 1))))
pgsim <- function(cl, R, n, a) {
    nw <- length(cl)
    RR <- rep(R / nw, nw)
    val <- do.call(rbind, parLapply(cl, RR, gsim, n, a))
    c(mn = mergeMeans(val[, "mn"], RR),
      md = mergeMeans(val[, "md"], RR),
      sd.mn = mergeSDs(val[, "sd.mn"], RR),
      sd.md = mergeSDs(val[, "sd.md"], RR))
}

## ## Load the parallel package, start a cluster of 2 worker processes,
## ## and initialize the random number streams
## > library(parallel)
## > cl <- makeCluster(2)
## > clusterSetRNGStream(cl, 123)

## ## Run some parallel simulations and compare to the sequential result.
## > system.time(val1 <- pgsim(cl, 100000, 10, 2))
##    user  system elapsed 
##   0.000   0.000   3.335 
## > val
##        mn        md     sd.mn     sd.md 
## 1.9981318 1.7298175 0.4474480 0.4827882 
## > val1
##        mn        md     sd.mn     sd.md 
## 1.9995171 1.7323369 0.4469518 0.4814748 
## > pgsim(cl, 100000, 10, 2)
##        mn        md     sd.mn     sd.md 
## 2.0010854 1.7328678 0.4469963 0.4815649

## ## reset the random number streams wit the same seed and check for
## ## reproducability
## > clusterSetRNGStream(cl, 123)
## > val2 <- pgsim(cl, 100000, 10, 2)
## > identical(val1, val2)

## ## shut down the cluster
## > stopCluster(cl)