pkg/R/plot.R

   1 # extractParam
   2 #
   3 # Extract successive values of a projection of the parameter(s)
   4 #
   5 # @inheritParams plotHist
   6 #
   7 extractParam <- function(mr, x=1, y=1)
   8 {
   9     # Obtain L vectors where L = number of res lists in mr
  10     lapply( mr, function(mr_list) {
  11         sapply(mr_list, function(m) m[x,y])
  12     } )
  13 }
  14
  15 #' plotHist
  16 #'
  17 #' Plot histogram
  18 #'
  19 #' @param mr Output of multiRun(), list of lists of functions results
  20 #' @param x Row index of the element inside the aggregated parameter
  21 #' @param y Column index of the element inside the aggregated parameter
  22 #'
  23 #' @examples
  24 #' \dontrun{
  25 #' β <- matrix(c(1,-2,3,1),ncol=2)
  26 #' mr <- multiRun(...) #see bootstrap example in ?multiRun : return lists of mu_hat
  27 #' μ <- normalize(β)
  28 #' for (i in 1:2)
  29 #'   mr[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
  30 #' plotHist(mr, 2, 1) #second row, first column}
  31 #' @export
  32 plotHist <- function(mr, x, y)
  33 {
  34     params <- extractParam(mr, x, y)
  35     L = length(params)
  36     # Plot histograms side by side
  37     par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
  38     for (i in 1:L)
  39         hist(params[[i]], breaks=40, freq=FALSE, xlab="Parameter value", ylab="Density")
  40 }
  41
  42 #' plotBox
  43 #'
  44 #' Draw boxplot
  45 #'
  46 #' @inheritParams plotHist
  47 #'
  48 #' @examples
  49 #' #See example in ?plotHist
  50 #' @export
  51 plotBox <- function(mr, x, y)
  52 {
  53     params <- extractParam(mr, x, y)
  54     L = length(params)
  55     # Plot boxplots side by side
  56     par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
  57     for (i in 1:L)
  58         boxplot(params[[i]], ylab="Parameter value")
  59 }
  60
  61 #' plotCoefs
  62 #'
  63 #' Draw coefs estimations + standard deviations
  64 #'
  65 #' @inheritParams plotHist
  66 #' @param params True value of parameters matrix
  67 #'
  68 #' @examples
  69 #' #See example in ?plotHist
  70 #' @export
  71 plotCoefs <- function(mr, params)
  72 {
  73     nf <- length(mr)
  74     L <- nrow(mr[[1]][[1]])
  75     K <- ncol(mr[[1]][[1]])
  76
  77     params_hat <- vector("list", nf)
  78     stdev <- vector("list", nf)
  79     for (i in 1:nf)
  80     {
  81         params_hat[[i]] <- matrix(nrow=L, ncol=K)
  82         stdev[[i]] <- matrix(nrow=L, ncol=K)
  83     }
  84     for (x in 1:L)
  85     {
  86         for (y in 1:K)
  87         {
  88             estims <- extractParam(mr, x, y)
  89             for (i in 1:nf)
  90             {
  91                 params_hat[[i]][x,y] <- mean(estims[[i]])
  92 #               stdev[[i]][x,y] <- sqrt( mean( (estims[[i]] - params[x,y])^2 ) )
  93                 # HACK remove extreme quantile in estims[[i]] before computing sd()
  94                 stdev[[i]][x,y] <- sd( estims[[i]] [ estims[[i]] < max(estims[[i]]) & estims[[i]] > min(estims[[i]]) ] )
  95             }
  96         }
  97     }
  98
  99     par(mfrow=c(1,nf), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
 100     params <- as.double(params)
 101     o <- order(params)
 102     for (i in 1:nf)
 103     {
 104         avg_param <- as.double(params_hat[[i]])
 105         std_param <- as.double(stdev[[i]])
 106         matplot(cbind(params[o],avg_param[o],avg_param[o]+std_param[o],avg_param[o]-std_param[o]),
 107             col=c(2,1,1,1), lty=c(1,1,2,2), type="l", lwd=2, xlab="param", ylab="value")
 108     }
 109
 110     #print(o) #not returning o to avoid weird Jupyter issue... (TODO:)
 111 }
 112
 113 #' plotQn
 114 #'
 115 #' Draw 3D map of objective function values
 116 #'
 117 #' @param N Number of starting points
 118 #' @param β Regression matrix (target)
 119 #' @param link Link function (logit or probit)
 120 #'
 121 #' @export
 122 plotQn <- function(N, n, p, β, b, link)
 123 {
 124     d <- nrow(β)
 125     K <- ncol(β)
 126     io <- generateSampleIO(n, p, β, b, link)
 127     op <- optimParams(K, link, list(X=io$X, Y=io$Y))
 128     # N random starting points gaussian (TODO: around true β?)
 129     res <- matrix(nrow=d*K+1, ncol=N)
 130     for (i in seq_len(N))
 131     {
 132         β_init <- rnorm(d*K)
 133         par <- op$run( c(rep(1/K,K-1), β_init, rep(0,K)) )
 134         par <- op$linArgs(par)
 135         Qn <- op$f(par)
 136         res[,i] = c(Qn, par[K:(K+d*K-1)])
 137     }
 138     res #TODO: plot this, not just return it...
 139 }