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