r_benchmark.R

system.time({
library(gtools)
library(MASS)
### data
yg <- galaxies
yg[78] <-26960
yg <- yg/1000
set.seed(123)
nsim <- 1e6
theta1 <- rnorm(nsim, 20, 10)
theta2 <- rnorm(nsim, 20, 10)
theta3 <- rnorm(nsim, 20, 10)
tau2 <- rgamma(nsim, 3, scale=1/20)
p <- rdirichlet(nsim, alpha=c(1,1,1))
l <- rep(1, nsim) # vector of likelihoods
for(y in yg) {
  l <- l* (p[,1]*dnorm(y, mean=theta1, sd=sqrt(1/tau2)) +
            p[,2]*dnorm(y, mean=theta2, sd=sqrt(1/tau2)) + 
            p[,3]*dnorm(y, mean=theta3, sd=sqrt(1/tau2)))
}
})

# user  system  elapsed

# work computer
# 36.28 3.03  39.42

# work linux server
# 33.523  0.177 34.741

# 15in rMBP mid 2015 (i7 2.2GHz)
# 15.065  1.253 16.338

# 13in MBA mid 2011 (i5 1.7GHz)
# 25.623  1.688 27.709

# 13in rMBP late 2013 (i5 2.4 GHz)
# 17.861  1.394 19.304

# nathan's hackintosh (i5 3.4 GHz)
# 14.461  0.982 15.440

# rpi 3
# 354.110 11.470  365.671

# hopkins biostats cluster (e5 1.7 GHz)
# 53.553  0.222 3.887