#
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])
- } )
+ # Obtain L vectors where L = number of res lists in mr
+ lapply( mr, function(mr_list) {
+ sapply(mr_list, function(m) m[x,y])
+ } )
}
#' plotHist
#' @export
plotHist <- function(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))
- for (i in 1:L)
- hist(params[[i]], breaks=40, freq=FALSE, xlab="Parameter value", ylab="Density")
+ 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))
+ for (i in 1:L)
+ hist(params[[i]], breaks=40, freq=FALSE, xlab="Parameter value", ylab="Density")
}
#' plotBox
#' @examples
#' #See example in ?plotHist
#' @export
-plotBox <- function(mr, x, y)
+plotBox <- function(mr, x, y, xtitle="")
{
- 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))
- for (i in 1:L)
- boxplot(params[[i]], ylab="Parameter value")
+ 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))
+ for (i in 1:L)
+ boxplot(params[[i]], xlab=xtitle, ylab="Parameter value")
}
#' plotCoefs
#'
#' @inheritParams plotHist
#' @param params True value of parameters matrix
+#' @param idx List index to process in mr
#'
#' @examples
#' #See example in ?plotHist
#' @export
-plotCoefs <- function(mr, params)
+plotCoefs <- function(mr, params, idx, xtitle="Parameter")
{
- nf <- length(mr)
- L <- nrow(mr[[1]][[1]])
- K <- ncol(mr[[1]][[1]])
+ L <- nrow(mr[[1]][[1]])
+ K <- ncol(mr[[1]][[1]])
- params_hat <- vector("list", nf)
- stdev <- vector("list", nf)
- for (i in 1:nf)
- {
- params_hat[[i]] <- matrix(nrow=L, ncol=K)
- stdev[[i]] <- matrix(nrow=L, ncol=K)
- }
- for (x in 1:L)
- {
- for (y in 1:K)
- {
- estims <- extractParam(mr, x, y)
- for (i in 1:nf)
- {
- params_hat[[i]][x,y] <- mean(estims[[i]])
-# stdev[[i]][x,y] <- sqrt( mean( (estims[[i]] - params[x,y])^2 ) )
- # HACK remove extreme quantile in estims[[i]] before computing sd()
- stdev[[i]][x,y] <- sd( estims[[i]] [ estims[[i]] < max(estims[[i]]) & estims[[i]] > min(estims[[i]]) ] )
- }
- }
- }
+ params_hat <- matrix(nrow=L, ncol=K)
+ stdev <- matrix(nrow=L, ncol=K)
+ for (x in 1:L)
+ {
+ 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 ) )
+ # 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]]) ] )
+ }
+ }
- par(mfrow=c(1,nf), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
- params <- as.double(params)
- o <- order(params)
- for (i in 1:nf)
- {
- avg_param <- as.double(params_hat[[i]])
- std_param <- as.double(stdev[[i]])
- 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,2,2), type="l", lwd=2, xlab="param", ylab="value")
- }
+ par(cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
+ params <- as.double(params)
+ 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="")
- #print(o) #not returning o to avoid weird Jupyter issue... (TODO:)
+ #print(o) #not returning o to avoid weird Jupyter issue... (TODO:)
}
#' plotQn
#' @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...
+ 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...
}