Fix multistart script
[morpheus.git] / reports / timings.R
CommitLineData
1d014a86
BA
1# flexmix optimization to get beta
2fmOptim <- 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)
15ourOptim <- 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
25getTimings <- 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
44link <- "logit"
45ncores <- 1
46N <- 100
47
48cmd_args <- commandArgs()
49for (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
64library(morpheus)
65library(flexmix)
66source("../patch_Bettina/FLXMRglm.R")
67
68tm <-
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 }
82tm_params <- list("link"=link, "N"=N, "nc"=ncores)
83
84save("tm", "tm_params", file="timings.RData")