| | 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") |
projects / morpheus.git / blame_incremental