pkg/R/plot.R

   1 # extractParam
   2 #
   3 # Extract successive values of a projection of the parameter(s).
   4 # The method works both on a list of lists of results,
   5 # or on a single list of parameters matrices.
   6 #
   7 # @inheritParams plotHist
   8 #
   9 .extractParam <- function(mr, x=1, y=1)
  10 {
  11   if (is.list(mr[[1]]))
  12   {
  13     # Obtain L vectors where L = number of res lists in mr
  14     return ( lapply( mr, function(mr_list) {
  15       sapply(mr_list, function(m) m[x,y])
  16     } ) )
  17   }
  18   sapply(mr, function(m) m[x,y])
  19 }
  20
  21 #' plotHist
  22 #'
  23 #' Plot compared histograms of a single parameter (scalar)
  24 #'
  25 #' @param mr Output of multiRun(), list of lists of functions results
  26 #' @param x Row index of the element inside the aggregated parameter
  27 #' @param y Column index of the element inside the aggregated parameter
  28 #' @param ... Additional graphical parameters (xlab, ylab, ...)
  29 #'
  30 #' @examples
  31 #' \donttest{
  32 #' β <- matrix(c(1,-2,3,1),ncol=2)
  33 #' mr <- multiRun(...) #see bootstrap example in ?multiRun
  34 #'                     #mr[[i]] is a list of estimated parameters matrices
  35 #' μ <- normalize(β)
  36 #' for (i in 1:2)
  37 #'   mr[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
  38 #' plotHist(mr, 2, 1) #second row, first column}
  39 #'
  40 #' @export
  41 plotHist <- function(mr, x, y, ...)
  42 {
  43   params <- .extractParam(mr, x, y)
  44   L = length(params)
  45   # Plot histograms side by side
  46   par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
  47   args <- list(...)
  48   for (i in 1:L)
  49   {
  50     hist(params[[i]], breaks=40, freq=FALSE,
  51       xlab=ifelse("xlab" %in% names(args), args$xlab, "Parameter value"),
  52       ylab=ifelse("ylab" %in% names(args), args$ylab, "Density"))
  53   }
  54 }
  55
  56 #' plotBox
  57 #'
  58 #' Draw compared boxplots of a single parameter (scalar)
  59 #'
  60 #' @inheritParams plotHist
  61 #'
  62 #' @examples
  63 #' \donttest{
  64 #' β <- matrix(c(1,-2,3,1),ncol=2)
  65 #' mr <- multiRun(...) #see bootstrap example in ?multiRun
  66 #'                     #mr[[i]] is a list of estimated parameters matrices
  67 #' μ <- normalize(β)
  68 #' for (i in 1:2)
  69 #'   mr[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
  70 #' plotBox(mr, 2, 1) #second row, first column}
  71 #'
  72 #' @export
  73 plotBox <- function(mr, x, y, ...)
  74 {
  75   params <- .extractParam(mr, x, y)
  76   L = length(params)
  77   # Plot boxplots side by side
  78   par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
  79   args <- list(...)
  80   for (i in 1:L)
  81   {
  82     boxplot(params[[i]],
  83       ifelse("ylab" %in% names(args), args$ylab, "Parameter value"))
  84   }
  85 }
  86
  87 #' plotCoefs
  88 #'
  89 #' Draw a graph of (averaged) coefficients estimations with their standard,
  90 #' deviations ordered by mean values.
  91 #' Note that the drawing does not correspond to a function; it is just a
  92 #' convenient way to visualize the estimated parameters.
  93 #'
  94 #' @param mr List of parameters matrices
  95 #' @param params True value of the parameters matrix
  96 #' @param ... Additional graphical parameters
  97 #'
  98 #' @examples
  99 #' \donttest{
 100 #' β <- matrix(c(1,-2,3,1),ncol=2)
 101 #' mr <- multiRun(...) #see bootstrap example in ?multiRun
 102 #'                     #mr[[i]] is a list of estimated parameters matrices
 103 #' μ <- normalize(β)
 104 #' for (i in 1:2)
 105 #'   mr[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
 106 #' params <- rbind( c(.5,.5), β, c(0,0) ) #p, β, b stacked in a matrix
 107 #' plotCoefs(mr[[1]], params)}
 108 #'
 109 #' @export
 110 plotCoefs <- function(mr, params, ...)
 111 {
 112   d <- nrow(mr[[1]])
 113   K <- ncol(mr[[1]])
 114
 115   params_hat <- matrix(nrow=d, ncol=K)
 116   stdev <- matrix(nrow=d, ncol=K)
 117   for (x in 1:d)
 118   {
 119     for (y in 1:K)
 120     {
 121       estims <- .extractParam(mr, x, y)
 122       params_hat[x,y] <- mean(estims)
 123       # Another way to compute stdev: using distances to true params
 124 #      stdev[x,y] <- sqrt( mean( (estims - params[x,y])^2 ) )
 125       # HACK remove extreme quantile in estims[[i]] before computing sd()
 126       stdev[x,y] <- sd(estims) #[ estims < max(estims) & estims > min(estims) ] )
 127     }
 128   }
 129
 130   par(cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
 131   params <- as.double(params)
 132   o <- order(params)
 133   avg_param <- as.double(params_hat)
 134   std_param <- as.double(stdev)
 135   args <- list(...)
 136   matplot(
 137     cbind(params[o],avg_param[o],
 138       avg_param[o]+std_param[o],avg_param[o]-std_param[o]),
 139     col=1, lty=c(1,5,3,3), type="l", lwd=2,
 140     xlab=ifelse("xlab" %in% names(args), args$xlab, "Parameter index"),
 141     ylab=ifelse("ylab" %in% names(args), args$ylab, "") )
 142
 143   #print(o) #not returning o to avoid weird Jupyter issue... (TODO:)
 144 }