X-Git-Url: https://git.auder.net/?p=morpheus.git;a=blobdiff_plain;f=pkg%2FR%2Fplot.R;h=da4fae53ac60504af51ca7b65be81a3df34c825b;hp=9727493e183cbdd4680d00c4a24c799a7120ce79;hb=2b3a6af5c55ac121405e3a8da721626ddf46b28b;hpb=2989133a738da99b78594cb00afdd3d87916b8bd

diff --git a/pkg/R/plot.R b/pkg/R/plot.R
index 9727493..da4fae5 100644
--- a/pkg/R/plot.R
+++ b/pkg/R/plot.R
@@ -1,90 +1,129 @@
 # extractParam
 #
-# Extract successive values of a projection of the parameter(s)
+# Extract successive values of a projection of the parameter(s).
+# The method works both on a list of lists of results,
+# or on a single list of parameters matrices.
 #
 # @inheritParams plotHist
 #
-extractParam <- function(mr, x=1, y=1)
+.extractParam <- function(mr, x=1, y=1)
 {
-  # Obtain L vectors where L = number of res lists in mr
-  lapply( mr, function(mr_list) {
-    sapply(mr_list, function(m) m[x,y])
-  } )
+  if (is.list(mr[[1]]))
+  {
+    # Obtain L vectors where L = number of res lists in mr
+    return ( lapply( mr, function(mr_list) {
+      sapply(mr_list, function(m) m[x,y])
+    } ) )
+  }
+  sapply(mr, function(m) m[x,y])
 }
 
 #' plotHist
 #'
-#' Plot histogram
+#' Plot compared histograms of a single parameter (scalar)
 #'
 #' @param mr Output of multiRun(), list of lists of functions results
 #' @param x Row index of the element inside the aggregated parameter
 #' @param y Column index of the element inside the aggregated parameter
+#' @param ... Additional graphical parameters (xlab, ylab, ...)
 #'
 #' @examples
 #' \donttest{
 #' β <- matrix(c(1,-2,3,1),ncol=2)
-#' mr <- multiRun(...) #see bootstrap example in ?multiRun : return lists of mu_hat
+#' mr <- multiRun(...) #see bootstrap example in ?multiRun
+#'                     #mr[[i]] is a list of estimated parameters matrices
 #' μ <- normalize(β)
 #' for (i in 1:2)
 #'   mr[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
 #' plotHist(mr, 2, 1) #second row, first column}
+#'
 #' @export
-plotHist <- function(mr, x, y)
+plotHist <- function(mr, x, y, ...)
 {
-  params <- extractParam(mr, x, y)
+  params <- .extractParam(mr, x, y)
   L = length(params)
   # Plot histograms side by side
   par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
+  args <- list(...)
   for (i in 1:L)
-    hist(params[[i]], breaks=40, freq=FALSE, xlab="Parameter value", ylab="Density")
+  {
+    hist(params[[i]], breaks=40, freq=FALSE,
+      xlab=ifelse("xlab" %in% names(args), args$xlab, "Parameter value"),
+      ylab=ifelse("ylab" %in% names(args), args$ylab, "Density"))
+  }
 }
 
 #' plotBox
 #'
-#' Draw boxplot
+#' Draw compared boxplots of a single parameter (scalar)
 #'
 #' @inheritParams plotHist
 #'
 #' @examples
-#' #See example in ?plotHist
+#' \donttest{
+#' β <- matrix(c(1,-2,3,1),ncol=2)
+#' mr <- multiRun(...) #see bootstrap example in ?multiRun
+#'                     #mr[[i]] is a list of estimated parameters matrices
+#' μ <- normalize(β)
+#' for (i in 1:2)
+#'   mr[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
+#' plotBox(mr, 2, 1) #second row, first column}
+#'
 #' @export
-plotBox <- function(mr, x, y, xtitle="")
+plotBox <- function(mr, x, y, ...)
 {
-  params <- extractParam(mr, x, y)
+  params <- .extractParam(mr, x, y)
   L = length(params)
   # Plot boxplots side by side
   par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
+  args <- list(...)
   for (i in 1:L)
-    boxplot(params[[i]], xlab=xtitle, ylab="Parameter value")
+  {
+    boxplot(params[[i]],
+      ifelse("ylab" %in% names(args), args$ylab, "Parameter value"))
+  }
 }
 
 #' plotCoefs
 #'
-#' Draw coefs estimations + standard deviations
+#' Draw a graph of (averaged) coefficients estimations with their standard,
+#' deviations ordered by mean values.
+#' Note that the drawing does not correspond to a function; it is just a
+#' convenient way to visualize the estimated parameters.
 #'
-#' @inheritParams plotHist
-#' @param params True value of parameters matrix
-#' @param idx List index to process in mr
+#' @param mr List of parameters matrices
+#' @param params True value of the parameters matrix
+#' @param ... Additional graphical parameters
 #'
 #' @examples
-#' #See example in ?plotHist
+#' \donttest{
+#' β <- matrix(c(1,-2,3,1),ncol=2)
+#' mr <- multiRun(...) #see bootstrap example in ?multiRun
+#'                     #mr[[i]] is a list of estimated parameters matrices
+#' μ <- normalize(β)
+#' for (i in 1:2)
+#'   mr[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
+#' params <- rbind( c(.5,.5), β, c(0,0) ) #p, β, b stacked in a matrix
+#' plotCoefs(mr[[1]], params)}
+#'
 #' @export
-plotCoefs <- function(mr, params, idx, xtitle="Parameter")
+plotCoefs <- function(mr, params, ...)
 {
-  L <- nrow(mr[[1]][[1]])
-  K <- ncol(mr[[1]][[1]])
+  d <- nrow(mr[[1]])
+  K <- ncol(mr[[1]])
 
-  params_hat <- matrix(nrow=L, ncol=K)
-  stdev <- matrix(nrow=L, ncol=K)
-  for (x in 1:L)
+  params_hat <- matrix(nrow=d, ncol=K)
+  stdev <- matrix(nrow=d, ncol=K)
+  for (x in 1:d)
   {
     for (y in 1:K)
     {
-      estims <- extractParam(mr, x, y)
-      params_hat[x,y] <- mean(estims[[idx]])
-#      stdev[x,y] <- sqrt( mean( (estims[[idx]] - params[x,y])^2 ) )
+      estims <- .extractParam(mr, x, y)
+      params_hat[x,y] <- mean(estims)
+      # Another way to compute stdev: using distances to true params
+#      stdev[x,y] <- sqrt( mean( (estims - params[x,y])^2 ) )
       # HACK remove extreme quantile in estims[[i]] before computing sd()
-      stdev[x,y] <- sd( estims[[idx]] ) #[ estims[[idx]] < max(estims[[idx]]) & estims[[idx]] > min(estims[[idx]]) ] )
+      stdev[x,y] <- sd(estims) #[ estims < max(estims) & estims > min(estims) ] )
     }
   }
 
@@ -93,39 +132,13 @@ plotCoefs <- function(mr, params, idx, xtitle="Parameter")
   o <- order(params)
   avg_param <- as.double(params_hat)
   std_param <- as.double(stdev)
-  matplot(cbind(params[o],avg_param[o],avg_param[o]+std_param[o],avg_param[o]-std_param[o]),
-    col=1, lty=c(1,5,3,3), type="l", lwd=2, xlab=xtitle, ylab="")
+  args <- list(...)
+  matplot(
+    cbind(params[o],avg_param[o],
+      avg_param[o]+std_param[o],avg_param[o]-std_param[o]),
+    col=1, lty=c(1,5,3,3), type="l", lwd=2,
+    xlab=ifelse("xlab" %in% names(args), args$xlab, "Parameter index"),
+    ylab=ifelse("ylab" %in% names(args), args$ylab, "") )
 
   #print(o) #not returning o to avoid weird Jupyter issue... (TODO:)
 }
-
-#' plotQn
-#'
-#' Draw 3D map of objective function values
-#'
-#' @param N Number of starting points
-#' @param n Number of points in sample
-#' @param p Vector of proportions
-#' @param b Vector of biases
-#' @param β Regression matrix (target)
-#' @param link Link function (logit or probit)
-#'
-#' @export
-plotQn <- function(N, n, p, β, b, link)
-{
-  d <- nrow(β)
-  K <- ncol(β)
-  io <- generateSampleIO(n, p, β, b, link)
-  op <- optimParams(K, link, list(X=io$X, Y=io$Y))
-  # N random starting points gaussian (TODO: around true β?)
-  res <- matrix(nrow=d*K+1, ncol=N)
-  for (i in seq_len(N))
-  {
-    β_init <- rnorm(d*K)
-    par <- op$run( c(rep(1/K,K-1), β_init, rep(0,K)) )
-    par <- op$linArgs(par)
-    Qn <- op$f(par)
-    res[,i] = c(Qn, par[K:(K+d*K-1)])
-  }
-  res #TODO: plot this, not just return it...
-}