reports/timings.R

   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 = list(p=rep(1/K,K-1), beta=mu, b=rep(0,K))
  20     optimParams(X, Y, K, link, M, 1)$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")