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