| 1 | # flexmix optimization to get beta |
| 2 | fmOptim <- function(X, Y, K, link) |
| 3 | { |
| 4 | dat <- as.data.frame( cbind(Y,X) ) |
| 5 | fm <- flexmix( cbind(Y, 1-Y) ~ .-Y, data=dat, k=K, |
| 6 | model = FLXMRglm(family = binomial(link = link)) ) |
| 7 | p <- mean(fm@posterior[["scaled"]][,1]) |
| 8 | out <- refit(fm) |
| 9 | beta_b <- sapply( seq_len(K), function(i) as.double( out@components[[1]][[i]][,1] ) ) |
| 10 | list("p"=p, "beta"=beta_b[2:nrow(beta_b),], "b"=beta_b[1,]) |
| 11 | NULL |
| 12 | } |
| 13 | |
| 14 | # Our package optimization for beta (using mu as a starting point) |
| 15 | ourOptim <- function(X, Y, K, link) |
| 16 | { |
| 17 | M <- computeMoments(X, Y) |
| 18 | mu <- computeMu(X, Y, list(K=K,M=M)) |
| 19 | x_init = c(1/2, as.double(mu), c(0,0)) |
| 20 | optimParams(K, link, list(M=M))$run(x_init) |
| 21 | NULL |
| 22 | } |
| 23 | |
| 24 | # Get timings for both methods with the same beta matrix |
| 25 | getTimings <- function(link) |
| 26 | { |
| 27 | timings <- list('fm'=matrix(0,nrow=10,ncol=7),'our'=matrix(0,nrow=10,ncol=7)) |
| 28 | K <- 2 |
| 29 | for (d in c(2,5,10)) |
| 30 | { |
| 31 | beta <- matrix(runif(d*K,min=-5,max=5),ncol=K) |
| 32 | for (logn in 4:6) |
| 33 | { |
| 34 | n <- 10^logn |
| 35 | io <- generateSampleIO(n, rep(1/K,K-1), beta, runif(K), link) |
| 36 | timings[['fm']][d,logn] <- system.time(fmOptim(io$X,io$Y,K,link))[3] |
| 37 | timings[['our']][d,logn] <- system.time(ourOptim(io$X,io$Y,K,link))[3] |
| 38 | } |
| 39 | } |
| 40 | timings |
| 41 | } |
| 42 | |
| 43 | #model = binomial |
| 44 | link <- "logit" |
| 45 | ncores <- 1 |
| 46 | N <- 100 |
| 47 | |
| 48 | cmd_args <- commandArgs() |
| 49 | for (arg in cmd_args) |
| 50 | { |
| 51 | if (substr(arg,1,1)!='-') |
| 52 | { |
| 53 | spl <- strsplit(arg,'=')[[1]] |
| 54 | if (spl[1] == "link") { |
| 55 | link <- spl[2] |
| 56 | } else if (spl[1] == "nc") { |
| 57 | ncores <- as.integer(spl[2]) |
| 58 | } else if (spl[1] == "N") { |
| 59 | N <- as.integer(spl[2]) |
| 60 | } |
| 61 | } |
| 62 | } |
| 63 | |
| 64 | library(morpheus) |
| 65 | library(flexmix) |
| 66 | source("../patch_Bettina/FLXMRglm.R") |
| 67 | |
| 68 | tm <- |
| 69 | if (ncores == 1) { |
| 70 | lapply(1:N, function(i) { |
| 71 | print(paste("Run",i)) |
| 72 | getTimings(link) |
| 73 | }) |
| 74 | } else { |
| 75 | library(parallel) |
| 76 | mclapply(1:N, function(i) { |
| 77 | print(paste("Run",i)) |
| 78 | getTimings(link) |
| 79 | }, |
| 80 | mc.preschedule=FALSE, mc.cores=ncores) |
| 81 | } |
| 82 | tm_params <- list("link"=link, "N"=N, "nc"=ncores) |
| 83 | |
| 84 | save("tm", "tm_params", file="timings.RData") |