From f4e42a2bc86f5b36a549f356033e4da8d07d0f81 Mon Sep 17 00:00:00 2001 From: Benjamin Auder Date: Mon, 16 Dec 2019 19:46:58 +0100 Subject: [PATCH] Reintroduce optional arg Mhat --- pkg/R/optimParams.R | 28 ++++--- reports/accuracy.R | 173 ++++++++++++++++++++++---------------------- 2 files changed, 104 insertions(+), 97 deletions(-) diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R index c1d7fe8..d8e2cf9 100644 --- a/pkg/R/optimParams.R +++ b/pkg/R/optimParams.R @@ -31,7 +31,7 @@ #' o$f( o$linArgs(par0) ) #' o$f( o$linArgs(par1) ) #' @export -optimParams <- function(X, Y, K, link=c("logit","probit")) +optimParams <- function(X, Y, K, link=c("logit","probit"), M=NULL) { # Check arguments if (!is.matrix(X) || any(is.na(X))) @@ -42,9 +42,19 @@ optimParams <- function(X, Y, K, link=c("logit","probit")) if (!is.numeric(K) || K!=floor(K) || K < 2) stop("K: integer >= 2") + if (is.null(M)) + { + # Precompute empirical moments + Mtmp <- computeMoments(X, Y) + M1 <- as.double(Mtmp[[1]]) + M2 <- as.double(Mtmp[[2]]) + M3 <- as.double(Mtmp[[3]]) + M <- c(M1, M2, M3) + } + # Build and return optimization algorithm object methods::new("OptimParams", "li"=link, "X"=X, - "Y"=as.integer(Y), "K"=as.integer(K)) + "Y"=as.integer(Y), "K"=as.integer(K), "Mhat"=as.double(M)) } #' Encapsulated optimization for p (proportions), β and b (regression parameters) @@ -82,18 +92,14 @@ setRefClass( "Check args and initialize K, d, W" callSuper(...) - if (!hasArg("X") || !hasArg("Y") || !hasArg("K") || !hasArg("li")) + if (!hasArg("X") || !hasArg("Y") || !hasArg("K") + || !hasArg("li") || !hasArg("Mhat")) + { stop("Missing arguments") - - # Precompute empirical moments - M <- computeMoments(X, Y) - M1 <- as.double(M[[1]]) - M2 <- as.double(M[[2]]) - M3 <- as.double(M[[3]]) - Mhat <<- c(M1, M2, M3) + } n <<- nrow(X) - d <<- length(M1) + d <<- ncol(X) W <<- diag(d+d^2+d^3) #initialize at W = Identity }, diff --git a/reports/accuracy.R b/reports/accuracy.R index ee08078..fd22a31 100644 --- a/reports/accuracy.R +++ b/reports/accuracy.R @@ -1,62 +1,63 @@ optimBeta <- function(N, n, K, p, beta, b, link, ncores) { - library(morpheus) - res <- multiRun( - list(n=n, p=p, beta=beta, b=b, K=K, link=link), - list( - # morpheus - function(fargs) { - library(morpheus) - K <- fargs$K - mu <- computeMu(fargs$X, fargs$Y, list(K=K)) - op <- optimParams(fargs$X, fargs$Y, K, fargs$link) + library(morpheus) + res <- multiRun( + list(n=n, p=p, beta=beta, b=b, K=K, link=link), + list( + # morpheus + function(fargs) { + library(morpheus) + K <- fargs$K + M <- computeMoments(fargs$X, fargs$Y) + mu <- computeMu(fargs$X, fargs$Y, list(K=K, M=M)) + op <- optimParams(fargs$X, fargs$Y, K, fargs$link, M) x_init <- list(p=rep(1/K,K-1), beta=mu, b=rep(0,K)) - res2 <- NULL - tryCatch({ + res2 <- NULL + tryCatch({ res2 <- do.call(rbind, op$run(x_init)) - }, error = function(e) {}) - res2 - } -# , -# # flexmix -# function(fargs) { -# library(flexmix) -# source("../patch_Bettina/FLXMRglm.R") -# K <- fargs$K -# dat <- as.data.frame( cbind(fargs$Y,fargs$X) ) -# res2 <- NULL -# tryCatch({ -# fm <- flexmix( cbind(V1, 1-V1) ~ .-V1, 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] ) -# } ) -# res2 <- rbind(p, beta_b[2:nrow(beta_b),], beta_b[1,]) -# }, error = function(e) { -# res2 <- NA -# }) -# res2 -# } - ), - prepareArgs = function(fargs, index) { - library(morpheus) - io = generateSampleIO(fargs$n, fargs$p, fargs$beta, fargs$b, fargs$link) - fargs$X = io$X - fargs$Y = io$Y - fargs - }, N=N, ncores=ncores, verbose=TRUE) - p <- c(p, 1-sum(p)) - for (i in 1:length(res)) { - for (j in N:1) { - if (is.null(res[[i]][[j]]) || is.na(res[[i]][[j]])) - res[[i]][[j]] <- NULL - } - print(paste("Count valid runs for ",i," = ",length(res[[i]]),sep="")) - res[[i]] <- alignMatrices(res[[i]], ref=rbind(p,beta,b), ls_mode="exact") - } - res + }, error = function(e) {}) + res2 + } +# , +# # flexmix +# function(fargs) { +# library(flexmix) +# source("../patch_Bettina/FLXMRglm.R") +# K <- fargs$K +# dat <- as.data.frame( cbind(fargs$Y,fargs$X) ) +# res2 <- NULL +# tryCatch({ +# fm <- flexmix( cbind(V1, 1-V1) ~ .-V1, 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] ) +# } ) +# res2 <- rbind(p, beta_b[2:nrow(beta_b),], beta_b[1,]) +# }, error = function(e) { +# res2 <- NA +# }) +# res2 +# } + ), + prepareArgs = function(fargs, index) { + library(morpheus) + io = generateSampleIO(fargs$n, fargs$p, fargs$beta, fargs$b, fargs$link) + fargs$X = io$X + fargs$Y = io$Y + fargs + }, N=N, ncores=ncores, verbose=TRUE) + p <- c(p, 1-sum(p)) + for (i in 1:length(res)) { + for (j in N:1) { + if (is.null(res[[i]][[j]]) || is.na(res[[i]][[j]])) + res[[i]][[j]] <- NULL + } + print(paste("Count valid runs for ",i," = ",length(res[[i]]),sep="")) + res[[i]] <- alignMatrices(res[[i]], ref=rbind(p,beta,b), ls_mode="exact") + } + res } #model = binomial; default values: @@ -69,46 +70,46 @@ ncores <- 1 cmd_args <- commandArgs() for (arg in cmd_args) { - if (substr(arg,1,1)!='-') { - spl <- strsplit(arg,'=')[[1]] - if (spl[1] == "nc") { - ncores <- as.integer(spl[2]) - } else if (spl[1] == "N") { - N <- as.integer(spl[2]) - } else if (spl[1] == "n") { - n <- as.integer(spl[2]) - } else if (spl[1] == "d") { - d <- as.integer(spl[2]) - } else if (spl[1] == "link") { - link <- spl[2] - } - } + if (substr(arg,1,1)!='-') { + spl <- strsplit(arg,'=')[[1]] + if (spl[1] == "nc") { + ncores <- as.integer(spl[2]) + } else if (spl[1] == "N") { + N <- as.integer(spl[2]) + } else if (spl[1] == "n") { + n <- as.integer(spl[2]) + } else if (spl[1] == "d") { + d <- as.integer(spl[2]) + } else if (spl[1] == "link") { + link <- spl[2] + } + } } if (d == 2) { - K <- 2 - p <- .5 - b <- c(-.2, .5) - beta <- matrix( c(1,-2, 3,1), ncol=K ) + K <- 2 + p <- .5 + b <- c(-.2, .5) + beta <- matrix( c(1,-2, 3,1), ncol=K ) } else if (d == 5) { - K <- 2 - p <- .5 - b <- c(-.2, .5) - beta <- matrix( c(1,2,-1,0,3, 2,-3,0,1,0), ncol=K ) + K <- 2 + p <- .5 + b <- c(-.2, .5) + beta <- matrix( c(1,2,-1,0,3, 2,-3,0,1,0), ncol=K ) } else if (d == 10) { - K <- 3 - p <- c(.3, .3) - b <- c(-.2, 0, .5) - beta <- matrix( c(1,2,-1,0,3,4,-1,-3,0,2, 2,-3,0,1,0,-1,-4,3,2,0, -1,1,3,-1,0,0,2,0,1,-2), ncol=K ) + K <- 3 + p <- c(.3, .3) + b <- c(-.2, 0, .5) + beta <- matrix( c(1,2,-1,0,3,4,-1,-3,0,2, 2,-3,0,1,0,-1,-4,3,2,0, -1,1,3,-1,0,0,2,0,1,-2), ncol=K ) } else if (d == 20) { - K <- 3 - p <- c(.3, .3) - b <- c(-.2, 0, .5) - beta <- matrix( c(1,2,-1,0,3,4,-1,-3,0,2,2,-3,0,1,0,-1,-4,3,2,0, -1,1,3,-1,0,0,2,0,1,-2,1,2,-1,0,3,4,-1,-3,0,2, 2,-3,0,1,0,-1,-4,3,2,0,1,1,2,2,-2,-2,3,1,0,0), ncol=K ) + K <- 3 + p <- c(.3, .3) + b <- c(-.2, 0, .5) + beta <- matrix( c(1,2,-1,0,3,4,-1,-3,0,2,2,-3,0,1,0,-1,-4,3,2,0, -1,1,3,-1,0,0,2,0,1,-2,1,2,-1,0,3,4,-1,-3,0,2, 2,-3,0,1,0,-1,-4,3,2,0,1,1,2,2,-2,-2,3,1,0,0), ncol=K ) } mr <- optimBeta(N, n, K, p, beta, b, link, ncores) mr_params <- list("N"=N, "nc"=ncores, "n"=n, "K"=K, "d"=d, "link"=link, - "p"=c(p,1-sum(p)), "beta"=beta, "b"=b) + "p"=c(p,1-sum(p)), "beta"=beta, "b"=b) save("mr", "mr_params", file=paste("res_",n,"_",d,"_",link,".RData",sep="")) -- 2.44.0