projects / morpheus.git / blame
| Commit | Line | Data |
|---|---|---|
| 1d014a86 BA |
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)) | |
| 5af71d43 BA |
19 | x_init = list(p=rep(1/K,K-1), beta=mu, b=rep(0,K)) |
| 20 | optimParams(X, Y, K, link, M, 1)$run(x_init) | |
| 1d014a86 BA |
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") |