Get index information into prepareArgs for multirun()
[morpheus.git] / 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 }