da4fae53ac60504af51ca7b65be81a3df34c825b
[morpheus.git] / 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 }