X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=reports%2Ftimings.R;fp=reports%2Ftimings.R;h=704fd82bd6c17695509b8f7f02ebef86a469e149;hb=476a79d4b332d378716e320fdfd95fa3e1d7be9f;hp=0000000000000000000000000000000000000000;hpb=0527116e18cd2ab53091de5c82d0d1d7356af4dc;p=morpheus.git

diff --git a/reports/timings.R b/reports/timings.R
new file mode 100644
index 0000000..704fd82
--- /dev/null
+++ b/reports/timings.R
@@ -0,0 +1,84 @@
+# flexmix optimization to get beta
+fmOptim <- function(X, Y, K, link)
+{
+	dat <- as.data.frame( cbind(Y,X) )
+	fm <- flexmix( cbind(Y, 1-Y) ~ .-Y, data=dat, k=K,
+		model = FLXMRglm(family = binomial(link = link)) )
+	p <- mean(fm@posterior[["scaled"]][,1])
+	out <- refit(fm)
+	beta_b <- sapply( seq_len(K), function(i) as.double( out@components[[1]][[i]][,1] ) )
+	list("p"=p, "beta"=beta_b[2:nrow(beta_b),], "b"=beta_b[1,])
+	NULL
+}
+
+# Our package optimization for beta (using mu as a starting point)
+ourOptim <- function(X, Y, K, link)
+{
+	M <- computeMoments(X, Y)
+	mu <- computeMu(X, Y, list(K=K,M=M))
+	x_init = c(1/2, as.double(mu), c(0,0))
+	optimParams(K, link, list(M=M))$run(x_init)
+	NULL
+}
+
+# Get timings for both methods with the same beta matrix
+getTimings <- function(link)
+{
+	timings <- list('fm'=matrix(0,nrow=10,ncol=7),'our'=matrix(0,nrow=10,ncol=7))
+	K <- 2
+	for (d in c(2,5,10))
+	{
+		beta <- matrix(runif(d*K,min=-5,max=5),ncol=K)
+		for (logn in 4:6)
+		{
+			n <- 10^logn
+			io <- generateSampleIO(n, rep(1/K,K-1), beta, runif(K), link)
+			timings[['fm']][d,logn] <- system.time(fmOptim(io$X,io$Y,K,link))[3]
+			timings[['our']][d,logn] <- system.time(ourOptim(io$X,io$Y,K,link))[3]
+		}
+	}
+	timings
+}
+
+#model = binomial
+link <- "logit"
+ncores <- 1
+N <- 100
+
+cmd_args <- commandArgs()
+for (arg in cmd_args)
+{
+	if (substr(arg,1,1)!='-')
+	{
+		spl <- strsplit(arg,'=')[[1]]
+		if (spl[1] == "link") {
+			link <- spl[2]
+		} else if (spl[1] == "nc") {
+			ncores <- as.integer(spl[2])
+		} else if (spl[1] == "N") {
+			N <- as.integer(spl[2])
+		}
+	}
+}
+
+library(morpheus)
+library(flexmix)
+source("../patch_Bettina/FLXMRglm.R")
+
+tm <-
+	if (ncores == 1) {
+		lapply(1:N, function(i) {
+			print(paste("Run",i))
+			getTimings(link)
+		})
+	} else {
+		library(parallel)
+		mclapply(1:N, function(i) {
+			print(paste("Run",i))
+			getTimings(link)
+		},
+		mc.preschedule=FALSE, mc.cores=ncores)
+	}
+tm_params <- list("link"=link, "N"=N, "nc"=ncores)
+
+save("tm", "tm_params", file="timings.RData")