Skip to content

Instantly share code, notes, and snippets.

@jacobcvt12
Last active Mar 3, 2017
Embed
What would you like to do?
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment