reports/accuracy.R

   1 optimBeta <- function(N, n, K, p, beta, b, link, ncores)
   2 {
   3     library(morpheus)
   4     res <- multiRun(
   5         list(n=n, p=p, beta=beta, b=b, K=K, link=link),
   6         list(
   7             # morpheus
   8             function(fargs) {
   9                 library(morpheus)
  10                 K <- fargs$K
  11                 mu <- computeMu(fargs$X, fargs$Y, list(K=K))
  12         op <- optimParams(fargs$X, fargs$Y, K, fargs$link)
  13         x_init <- list(p=rep(1/K,K-1), beta=mu, b=rep(0,K))
  14                 res2 <- NULL
  15                 tryCatch({
  16           res2 <- do.call(rbind, op$run(x_init))
  17                 }, error = function(e) {})
  18                 res2
  19             }
  20 #           ,
  21 #           # flexmix
  22 #           function(fargs) {
  23 #               library(flexmix)
  24 #               source("../patch_Bettina/FLXMRglm.R")
  25 #               K <- fargs$K
  26 #               dat <- as.data.frame( cbind(fargs$Y,fargs$X) )
  27 #               res2 <- NULL
  28 #               tryCatch({
  29 #                   fm <- flexmix( cbind(V1, 1-V1) ~ .-V1, data=dat, k=K,
  30 #                       model = FLXMRglm(family = binomial(link = link)) )
  31 #                   p <- mean(fm@posterior[["scaled"]][,1])
  32 #                   out <- refit(fm)
  33 #                   beta_b <- sapply( seq_len(K), function(i) {
  34 #                       as.double( out@components[[1]][[i]][,1] )
  35 #                   } )
  36 #                   res2 <- rbind(p, beta_b[2:nrow(beta_b),], beta_b[1,])
  37 #               }, error = function(e) {
  38 #                   res2 <- NA
  39 #               })
  40 #               res2
  41 #           }
  42         ),
  43         prepareArgs = function(fargs, index) {
  44             library(morpheus)
  45             io = generateSampleIO(fargs$n, fargs$p, fargs$beta, fargs$b, fargs$link)
  46             fargs$X = io$X
  47             fargs$Y = io$Y
  48             fargs
  49         }, N=N, ncores=ncores, verbose=TRUE)
  50     p <- c(p, 1-sum(p))
  51     for (i in 1:length(res)) {
  52         for (j in N:1) {
  53             if (is.null(res[[i]][[j]]) || is.na(res[[i]][[j]]))
  54                 res[[i]][[j]] <- NULL
  55         }
  56         print(paste("Count valid runs for ",i," = ",length(res[[i]]),sep=""))
  57         res[[i]] <- alignMatrices(res[[i]], ref=rbind(p,beta,b), ls_mode="exact")
  58     }
  59     res
  60 }
  61
  62 #model = binomial; default values:
  63 link = "logit"
  64 N <- 10
  65 d <- 2
  66 n <- 1e4
  67 ncores <- 1
  68
  69 cmd_args <- commandArgs()
  70 for (arg in cmd_args)
  71 {
  72     if (substr(arg,1,1)!='-') {
  73         spl <- strsplit(arg,'=')[[1]]
  74         if (spl[1] == "nc") {
  75             ncores <- as.integer(spl[2])
  76         } else if (spl[1] == "N") {
  77             N <- as.integer(spl[2])
  78         } else if (spl[1] == "n") {
  79             n <- as.integer(spl[2])
  80         } else if (spl[1] == "d") {
  81             d <- as.integer(spl[2])
  82         } else if (spl[1] == "link") {
  83             link <- spl[2]
  84         }
  85     }
  86 }
  87
  88 if (d == 2) {
  89     K <- 2
  90     p <- .5
  91     b <- c(-.2, .5)
  92     beta <- matrix( c(1,-2, 3,1), ncol=K )
  93 } else if (d == 5) {
  94     K <- 2
  95     p <- .5
  96     b <- c(-.2, .5)
  97     beta <- matrix( c(1,2,-1,0,3, 2,-3,0,1,0), ncol=K )
  98 } else if (d == 10) {
  99     K <- 3
 100     p <- c(.3, .3)
 101     b <- c(-.2, 0, .5)
 102     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 )
 103 } else if (d == 20) {
 104     K <- 3
 105     p <- c(.3, .3)
 106     b <- c(-.2, 0, .5)
 107     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 )
 108 }
 109
 110 mr <- optimBeta(N, n, K, p, beta, b, link, ncores)
 111 mr_params <- list("N"=N, "nc"=ncores, "n"=n, "K"=K, "d"=d, "link"=link,
 112     "p"=c(p,1-sum(p)), "beta"=beta, "b"=b)
 113
 114 save("mr", "mr_params", file=paste("res_",n,"_",d,"_",link,".RData",sep=""))