reports/accuracy.R

   1 optimBeta <- function(N, n, K, p, beta, b, link, weights, ncores)
   2 {
   3     library(morpheus)
   4     res <- multiRun(
   5         list(n=n, p=p, beta=beta, b=b, optargs=list(K=K, link=link, weights=weights)),
   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 weights <- c(1,1,1)
  75
  76 cmd_args <- commandArgs()
  77 for (arg in cmd_args)
  78 {
  79     if (substr(arg,1,1)!='-') {
  80         spl <- strsplit(arg,'=')[[1]]
  81         if (spl[1] == "nc") {
  82             ncores <- as.integer(spl[2])
  83         } else if (spl[1] == "N") {
  84             N <- as.integer(spl[2])
  85         } else if (spl[1] == "n") {
  86             n <- as.integer(spl[2])
  87         } else if (spl[1] == "d") {
  88             d <- as.integer(spl[2])
  89         } else if (spl[1] == "link") {
  90             link <- spl[2]
  91         } else if (spl[1] == "weights") {
  92           weights <- unlist(strsplit(spl[2], ","))
  93         }
  94     }
  95 }
  96
  97 if (d == 2) {
  98     K <- 2
  99     p <- .5
 100     b <- c(-.2, .5)
 101     beta <- matrix( c(1,-2, 3,1), ncol=K )
 102 } else if (d == 5) {
 103     K <- 2
 104     p <- .5
 105     b <- c(-.2, .5)
 106     beta <- matrix( c(1,2,-1,0,3, 2,-3,0,1,0), ncol=K )
 107 } else if (d == 10) {
 108     K <- 3
 109     p <- c(.3, .3)
 110     b <- c(-.2, 0, .5)
 111     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 )
 112 } else if (d == 20) {
 113     K <- 3
 114     p <- c(.3, .3)
 115     b <- c(-.2, 0, .5)
 116     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 )
 117 }
 118
 119 mr <- optimBeta(N, n, K, p, beta, b, link, weights, ncores)
 120 mr_params <- list("N"=N, "n"=n, "K"=K, "d"=d, "link"=link,
 121     "p"=c(p,1-sum(p)), "beta"=beta, "b"=b, "weights"=weights)
 122
 123 save("mr", "mr_params", file=paste("multirun_",n,"_",d,"_",link,".RData",sep=""))