#' Elbow plot of the loss function.
#'
#' Boxplots of the output \code{L} from \code{\link{cve}} over \code{k} from
#' \code{min.dim} to \code{max.dim}. For given \code{k}, \code{L} corresponds
#' to \eqn{L_n(V, X_i)} where \eqn{V} is the minimizer of \eqn{L_n(V)} where
#' \eqn{V} is an element of a Stiefel manifold (see
#' Fertl, L. and Bura, E. (2019)).
#'
#' @param x an object of class \code{"cve"}, usually, a result of a call to
#'          \code{\link{cve}} or \code{\link{cve.call}}.
#' @param ... Pass through parameters to [\code{\link{plot}}] and
#'      [\code{\link{lines}}]
#'
#' @examples
#' # create B for simulation
#' B <- cbind(rep(1, 6), (-1)^seq(6)) / sqrt(6)
#' 
#' set.seed(21)
#' # creat predictor data x ~ N(0, I_p)
#' X <- matrix(rnorm(600), 100)
#' 
#' # simulate response variable
#' #    y = f(B'x) + err
#' # with f(x1, x2) = x1^2 + 2 x2 and err ~ N(0, 0.25^2)
#' Y <- (X %*% B[, 1])^2 + 2 * X %*% B[, 2] + rnorm(100, 0, .1)
#'
#' # Create bandwidth estimation function
#' estimate.bandwidth <- function(X, k, nObs) {
#'     n <- nrow(X)
#'     p <- ncol(X)
#'     X_c <- scale(X, center = TRUE, scale = FALSE)
#'     2 * qchisq((nObs - 1) / (n - 1), k) * sum(X_c^2) / (n * p)
#' }
#' # calculate cve with method 'simple' for k = min.dim,...,max.dim
#' cve.obj.simple <- cve(Y ~ X, h = estimate.bandwidth, nObs = sqrt(nrow(X)))
#' 
#' # elbow plot
#' plot(cve.obj.simple)
#'
#' @references Fertl, L. and Bura, E. (2019), Conditional Variance
#' Estimation for Sufficient Dimension Reduction. Working Paper.
#'
#' @seealso see \code{\link{par}} for graphical parameters to pass through
#'      as well as \code{\link{plot}}, the standard plot utility.
#' @method plot cve
#' @importFrom graphics plot lines points boxplot
#' @export
plot.cve <- function(x, ...) {
    L <- c()
    k <- c()
    for (dr.k in x$res) {
        if (class(dr.k) == 'cve.k') {
            k <- c(k, as.character(dr.k$k))
            L <- c(L, dr.k$L)
        }
    }
    L <- matrix(L, ncol = length(k)) / var(x$Y)
    boxplot(L, main = "elbow plot",
            xlab = "SDR dimension",
            ylab = "Sample loss distribution",
            names = k)
}