.packageName <- "rrcov"
##  rrcov : Scalable Robust Estimators with High Breakdown Point
##
##  This program is free software; you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation; either version 2 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program; if not, write to the Free Software
##  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
##
##  I would like to thank Peter Rousseeuw and Katrien van Driessen for 
##  providing the initial code of this function.

covMcd <- function(x, 
                   cor=FALSE, 
                   alpha=1/2, 
                   nsamp=500, 
                   seed=0, 
                   print.it=FALSE)
{
    quan.f <- function(alpha, n, rk)
    {
        quan <- floor(2 * floor((n+rk+1)/2) - n + 2 * (n - floor((n+rk+1)/2)) * alpha)
        return(quan)
    }

    correctiefactor.s <- function(p, n, alpha)
    {
        if(p > 2) {
            coeffqpkwad875 <- matrix(c(-0.455179464070565, 
                1.11192541278794, 2, -0.294241208320834, 
                1.09649329149811, 3), ncol = 2, byrow = FALSE)
            dimnames(coeffqpkwad875) <- list(c("alfaq", "betaq", 
                "qwaarden"), c("coeffqpkwad875.q2", 
                "coeffqpkwad875.q3"))
            coeffqpkwad500 <- matrix(c(-1.42764571687802, 
                1.26263336932151, 2, -1.06141115981725, 
                1.28907991440387, 3), ncol = 2, byrow = FALSE)
            dimnames(coeffqpkwad500) <- list(c("alfaq", "betaq", 
                "qwaarden"), c("coeffqpkwad500.q2", 
                "coeffqpkwad500.q3"))
            y1.500 <- 1 + (coeffqpkwad500[1, 1] * 1)/p^{
                coeffqpkwad500[2, 1]
            }
            y2.500 <- 1 + (coeffqpkwad500[1, 2] * 1)/p^{
                coeffqpkwad500[2, 2]
            }
            y1.875 <- 1 + (coeffqpkwad875[1, 1] * 1)/p^{
                coeffqpkwad875[2, 1]
            }
            y2.875 <- 1 + (coeffqpkwad875[1, 2] * 1)/p^{
                coeffqpkwad875[2, 2]
            }
            y1.500 <- log(1 - y1.500)
            y2.500 <- log(1 - y2.500)
            y.500 <- c(y1.500, y2.500)
            A.500 <- matrix(c(1, log(1/(coeffqpkwad500[3, 1] * p^2)
                ), 1, log(1/(coeffqpkwad500[3, 2] * p^2))), 
                ncol = 2, byrow = TRUE)
            coeffic.500 <- solve(A.500, y.500)
            y1.875 <- log(1 - y1.875)
            y2.875 <- log(1 - y2.875)
            y.875 <- c(y1.875, y2.875)
            A.875 <- matrix(c(1, log(1/(coeffqpkwad875[3, 1] * p^2)
                ), 1, log(1/(coeffqpkwad875[3, 2] * p^2))), 
                ncol = 2, byrow = TRUE)
            coeffic.875 <- solve(A.875, y.875)
            fp.500.n <- 1 - (exp(coeffic.500[1]) * 1)/n^{
                coeffic.500[2]
            }
            fp.875.n <- 1 - (exp(coeffic.875[1]) * 1)/n^{
                coeffic.875[2]
            }
        }
        else {
            if(p == 2) {
                fp.500.n <- 1 - (exp(0.673292623522027) * 1)/
                  n^{
                  0.691365864961895
                }
                fp.875.n <- 1 - (exp(0.446537815635445) * 1)/
                  n^{
                  1.06690782995919
                }
            }
            if(p == 1) {
                fp.500.n <- 1 - (exp(0.262024211897096) * 1)/
                  n^{
                  0.604756680630497
                }
                fp.875.n <- 1 - (exp(-0.351584646688712) * 1)/
                  n^{
                  1.01646567502486
                }
            }
        }
        if((0.5 <= alpha) && (alpha <= 0.875)) {
            fp.alpha.n <- fp.500.n + (fp.875.n - fp.500.n)/0.375 * (
                alpha - 0.5)
        }
        if((0.875 < alpha) && (alpha <= 1)) {
            fp.alpha.n <- fp.875.n + (1 - fp.875.n)/0.125 * (alpha - 
                0.875)
        }
        return(1/fp.alpha.n)
    }
    
correctiefactor.rew.s <- function(p, n, alpha)
    {
        if(p > 2) {
            coeffrewqpkwad875 <- matrix(c(-0.544482443573914 , 
                 1.25994483222292   , 2, -0.343791072183285 , 
                1.25159004257133 , 3), ncol = 2, byrow = FALSE)
            dimnames(coeffrewqpkwad875) <- list(c("alfaq", "betaq", 
                "qwaarden"), c("coeffrewqpkwad875.q2", 
                "coeffrewqpkwad875.q3"))
            coeffrewqpkwad500 <- matrix(c(-1.02842572724793 , 
                1.67659883081926 , 2, -0.26800273450853 , 
                1.35968562893582 , 3), ncol = 2, byrow = FALSE)
            dimnames(coeffrewqpkwad500) <- list(c("alfaq", "betaq", 
                "qwaarden"), c("coeffrewqpkwad500.q2", 
                "coeffrewqpkwad500.q3"))
            y1.500 <- 1 + (coeffrewqpkwad500[1, 1] * 1)/p^{
                coeffrewqpkwad500[2, 1]
            }
            y2.500 <- 1 + (coeffrewqpkwad500[1, 2] * 1)/p^{
                coeffrewqpkwad500[2, 2]
            }
            y1.875 <- 1 + (coeffrewqpkwad875[1, 1] * 1)/p^{
                coeffrewqpkwad875[2, 1]
            }
            y2.875 <- 1 + (coeffrewqpkwad875[1, 2] * 1)/p^{
                coeffrewqpkwad875[2, 2]
            }
            y1.500 <- log(1 - y1.500)
            y2.500 <- log(1 - y2.500)
            y.500 <- c(y1.500, y2.500)
            A.500 <- matrix(c(1, log(1/(coeffrewqpkwad500[3, 1] * p^
                2)), 1, log(1/(coeffrewqpkwad500[3, 2] * p^2))),
                ncol = 2, byrow = TRUE)
            coeffic.500 <- solve(A.500, y.500)
            y1.875 <- log(1 - y1.875)
            y2.875 <- log(1 - y2.875)
            y.875 <- c(y1.875, y2.875)
            A.875 <- matrix(c(1, log(1/(coeffrewqpkwad875[3, 1] * p^
                2)), 1, log(1/(coeffrewqpkwad875[3, 2] * p^2))),
                ncol = 2, byrow = TRUE)
            coeffic.875 <- solve(A.875, y.875)
            fp.500.n <- 1 - (exp(coeffic.500[1]) * 1)/n^{
                coeffic.500[2]
            }
            fp.875.n <- 1 - (exp(coeffic.875[1]) * 1)/n^{
                coeffic.875[2]
            }
        }
        else {
            if(p == 2) {
                fp.500.n <- 1 - (exp( 3.11101712909049  ) * 1)/n^
                  {
                  1.91401056721863 
                }
                fp.875.n <- 1 - (exp( 0.79473550581058  ) * 1)/
                  n^{
                  1.10081930350091 
                }
            }
            if(p == 1) {
                fp.500.n <- 1 - (exp( 1.11098143415027  ) * 1)/n^
                  {
                  1.5182890270453 
                }
                fp.875.n <- 1 - (exp( -0.66046776772861 ) * 1)/n^
                  {
                 0.88939595831888  
                }
            }
        }
        if((0.5 <= alpha) && (alpha <= 0.875)) {
            fp.alpha.n <- fp.500.n + (fp.875.n - fp.500.n)/0.375 * (
                alpha - 0.5)
        }
        if((0.875 < alpha) && (alpha <= 1)) {
            fp.alpha.n <- fp.875.n + (1 - fp.875.n)/0.125 * (alpha - 
                0.875)
        }
        return(1/fp.alpha.n)
    }

    # vt:: not necessary for R 1.8 and higher - the R function determinant() is used
    #
    # Determinant for square real matrix x. Adapted from det.Matrix in library(Matrix).
    #
    #    determinant <- function(x, logarithm = TRUE){
    #...

    # vt:: tolerance to be used for computing the mahalanobis distances (default = 1e-7)
    tol = 1e-10

    if(is.vector(x) || (is.matrix(x) && !is.data.frame(x))) {
        if(!is.numeric(x))
            stop(message = "x is not a numeric dataframe or matrix.")
    }
    if((!is.vector(x) && !is.matrix(x)) || is.data.frame(x)) {
        if((!is.data.frame(x) && !is.numeric(x)) || (!all(sapply(x,data.class) == "numeric")))
            stop(message = "x is not a numeric dataframe or matrix.")
    }
    
    #vt:: if the data is supplied as a data.frame, the following expressions results in an error
    # as workaround convert the data.frame to a matrix
    if(is.data.frame(x))
        x <- as.matrix(x)
    
    if(!is.matrix(x))
        x <- array(x, c(length(x), 1), list(names(x), deparse(
            substitute(data))))
    x <- as.matrix(x)
    dimn <- dimnames(x)
    na.x <- !is.finite(x %*% rep(1, ncol(x)))
    ok <- !na.x
    x <- x[ok,  , drop = FALSE]
    dx <- dim(x)
    if(!length(dx))
        stop("All observations have missing values!")
    n <- dx[1]
    p <- dx[2]
    if(n < 2 * p)
        stop("Need at least 2*(number of variables) observations ")
    jmin <- floor((n + p + 1)/2)
    if(alpha < 1/2) {
        stop("The MCD must cover at least", jmin, "observations")
    }
    else if(alpha > 1)
        stop("alpha is out of range")
    quan <- quan.f(alpha, n, p)

    # Compute the classical estimates - alpha=1
    if(alpha == 1) {
        mcd <- cov.wt(x)$cov
        loc <- as.vector(apply(x, 2, mean))
        obj <- determinant(mcd, log = TRUE)$modulus[1]
        if(( - obj/p) > 50) {
            ans <- list()
            ans$cov <- mcd
            dimnames(ans$cov) <- list(dimn[[2]], dimn[[2]])
            ans$center <- loc
            if(length(dimn[[2]]))
                names(ans$center) <- dimn[[2]]
            ans$n.obs <- n
            ans$call <- match.call()
            ans$method <- paste(
                "Minimum Covariance Determinant Estimator.")
            ans$method <- paste(ans$method, 
                "\nThe classical covariance matrix is singular."
                )
            if(!print.it) {
                cat("The classical covariance matrix is singular.\n"
                  )
            }
            ans$alpha <- alpha
            ans$quan <- quan
            ans$raw.cov <- mcd
            dimnames(ans$raw.cov) <- list(dimn[[2]], dimn[[2]])
            ans$raw.center <- loc
            if(length(dimn[[2]]))
                names(ans$raw.center) <- dimn[[2]]
            ans$crit <- exp(obj)
            ans$mcd.wt <- rep(NA, length(na.x))
            ans$mcd.wt[ok] <- rep(1, sum(ok == TRUE))
        }
        else {
            mah <- mahalanobis(x, loc, mcd, tol.inv=tol)        # VT:: 01.09.2004 - bug in alpha=1 
                                                                # (tol instead of tol.inv as parameter name)
            weights <- ifelse(mah < qchisq(0.975, p), 1, 0)
            ans <- cov.wt(x, wt = weights, cor)
            ans$cov <- sum(weights)/(sum(weights) - 1) * ans$cov    
        
            #Consistency factor for reweighted MCD
            if(sum(weights) == n)
                cdelta.rew <- 1
            else {
                qdelta.rew <- qchisq(sum(weights)/n, p)
                cdeltainvers.rew <- pgamma(qdelta.rew/2, p/2 + 
                  1)/(sum(weights)/n)
                cdelta.rew <- 1/cdeltainvers.rew
            }
            ans$cov <- ans$cov * cdelta.rew
            ans$call <- match.call()
            ans$method <- paste("Minimum Covariance Determinant Estimator.")
            if( - (determinant(ans$cov, log = TRUE)$modulus[1] - 0)/p > 50) {
                ans$method <- paste(ans$method, "\nThe reweighted MCD scatter matrix is singular.")
                if(!print.it) {
                  cat("The reweighted MCD scatter matrix is singular.\n")
                }
                ans$alpha <- alpha
                ans$quan <- quan
                ans$raw.cov <- mcd
                dimnames(ans$raw.cov) <- list(dimn[[2]], dimn[[2]])
                ans$raw.center <- loc
                if(length(dimn[[2]]))
                    names(ans$raw.center) <- dimn[[2]]
                ans$crit <- exp(obj)
                ans$mcd.wt <- rep(NA, length(na.x))
                ans$mcd.wt[ok] <- weights
                if(length(dimn[[1]]))
                    names(ans$mcd.wt) <- dimn[[1]]
                ans$wt <- NULL
                ans$X <- x
                if(length(dimn[[1]]))
                    dimnames(ans$X)[[1]] <- names(ans$mcd.wt)[ok]
                else {
                    xx <- seq(1, length(na.x))
                    dimnames(ans$X) <- list(NULL, NULL)
                    dimnames(ans$X)[[1]] <- xx[ok]
                }
                ans$method <- paste(ans$method, "\nThe minimum covariance determinant estimates based on", n, 
                    "observations \nare equal to the classical estimates.")
                if(print.it) {
                  cat(ans$method, "\n")
                }
                class(ans) <- "mcd"
                attr(ans, "call") <- sys.call()
                return(ans)
            }
            else {
                ans$alpha <- alpha
                ans$quan <- quan
                ans$raw.cov <- mcd
                dimnames(ans$raw.cov) <- list(dimn[[2]], dimn[[2]])
                ans$raw.center <- loc
                if(length(dimn[[2]]))
                    names(ans$raw.center) <- dimn[[2]]
                ans$crit <- exp(obj)
                mah <- mahalanobis(x, ans$center, ans$cov, tol.inv=tol)
            }
            ans$mcd.wt <- rep(NA, length(na.x))
            ans$mcd.wt[ok] <- ifelse(mah < qchisq(0.975, p), 1, 0)
        }
        if(length(dimn[[1]]))
            names(ans$mcd.wt) <- dimn[[1]]
        ans$wt <- NULL
        ans$X <- x
        if(length(dimn[[1]]))
            dimnames(ans$X)[[1]] <- names(ans$mcd.wt)[ok]
        else {
            xx <- seq(1, length(na.x))
            dimnames(ans$X) <- list(NULL, NULL)
            dimnames(ans$X)[[1]] <- xx[ok]
        }
        ans$method <- paste(ans$method, 
            "\nThe minimum covariance determinant estimates based on",
            n, "observations \nare equal to the classical estimates."
            )
        if(print.it) {
            cat(ans$method, "\n")
        }
        class(ans) <- "mcd"
        attr(ans, "call") <- sys.call()
        return(ans)
    }   #end alpha=1
    
    storage.mode(x) <- "double"
    storage.mode(quan) <- "integer"
    initcov <- matrix(0, nrow = p * p, ncol = 1)
    adcov <- matrix(0, nrow = p * p, ncol = 1)
    initmean <- matrix(0, nrow = p, ncol = 1)
    inbest <- matrix(10000, nrow = quan, ncol = 1)
    plane <- matrix(0, nrow = 5, ncol = p)
    deter <- 0
    weights <- matrix(0, nrow = n, ncol = 1)
    fit <- 0
    kount <- 0
    storage.mode(n) <- "integer"
    storage.mode(p) <- "integer"
    storage.mode(nsamp) <- "integer"
    storage.mode(initcov) <- "double"
    storage.mode(adcov) <- "double"
    storage.mode(initmean) <- "double"
    storage.mode(inbest) <- "integer"
    storage.mode(plane) <- "double"
    storage.mode(deter) <- "double"
    storage.mode(weights) <- "integer"
    storage.mode(fit) <- "integer"
    storage.mode(kount) <- "integer"
    storage.mode(seed) <- "integer"
    mcd <- .Fortran("rffastmcd",
        x,
        n,
        p,
        quan,
        nsamp,
        initcovariance = initcov,
        initmean = initmean,
        best=inbest,
        mcdestimate = deter,
        weights = weights,
        exactfit = fit,
        coeff = plane,
        kount = kount,
        adjustcov = adcov,
        seed,
        PACKAGE="rrcov")  
    
    # Compute the consistency correction factor for the raw MCD (see calfa in croux and haesbroeck)
    qalpha <- qchisq(quan/n, p)
    calphainvers <- pgamma(qalpha/2, p/2 + 1)/(quan/n)
    calpha <- 1/calphainvers
    correct <- correctiefactor.s(p, n, alpha)
    
    if(p == 1) 
    {
        scale <- sqrt(calpha) * as.double(mcd$initcovariance) * sqrt(correct)
    }else 
    {
        # Apply correction factor to the raw estimates and use them to compute weights
        mcd$initcovariance <- calpha * mcd$initcovariance * correct
        dim(mcd$initcovariance) <- c(p, p)
        if(mcd$exactfit == 0) 
        {
            mah <- mahalanobis(x, mcd$initmean, mcd$initcovariance, tol.inv = tol)
            mcd$weights <- ifelse(mah < qchisq(0.975, p), 1, 0)
        }
    }
    
    #The number of variables is 1 - compute univariate location and scale estimates
    if(p == 1) 
    {
        center <- as.double(mcd$initmean)
        if(abs(scale - 0) < 1e-07) 
        {
            ans <- list()
            ans$cov <- 0
            names(ans$cov) <- dimn[[2]][1]
            ans$center <- center
            names(ans$center) <- dimn[[2]][1]
            ans$n.obs <- n
            ans$call <- match.call()    
            ans$method <- paste(
                "Univariate location and scale estimation.\nMore than",
                quan, "of the observations are identical.")
            ans$alpha <- alpha
            ans$quan <- quan
            ans$raw.cov <- 0
            names(ans$raw.cov) <- dimn[[2]][1]
            ans$raw.center <- center
            names(ans$raw.center) <- dimn[[2]][1]
            ans$crit <- 0  
            ans$mcd.wt <- rep(NA, length(na.x))
            ans$mcd.wt[ok] <- as.vector(ifelse(abs(x - center) <    1e-07, 1, 0))
            if(length(dimn[[1]]))
                names(ans$mcd.wt) <- dimn[[1]]
            if(print.it) {
                cat(ans$method, "\n")
            }
            ans$X <- x
            if(length(dimn[[1]]))
                dimnames(ans$X)[[1]] <- names(ans$mcd.wt)[ok]
            else {
                xx <- seq(1, length(na.x))
                dimnames(ans$X) <- list(NULL, NULL)
                dimnames(ans$X)[[1]] <- xx[ok]
            }
            
            class(ans) <- "mcd"
            attr(ans, "call") <- sys.call()
            return(ans)
        }
        
        # Compute the weights for the raw MCD in case p=1
        quantiel <- qchisq(0.975,p)
        weights <- ifelse(((x - center)/scale)^2  < quantiel, 1, 0)
        ans <- cov.wt(x, wt = weights, cor = cor)
        ans$cov <- sum(weights)/(sum(weights) - 1) * ans$cov
        
        #Apply the correction factor for the reweighted cov
        if(sum(weights) == n)
        {
            cdelta.rew <- 1
            correct.rew <- 1
        }else 
        {
            qdelta.rew <- qchisq(sum(weights)/n, p)
            cdeltainvers.rew <- pgamma(qdelta.rew/2, p/2 + 1)/(sum(weights)/n)
            cdelta.rew <- 1/cdeltainvers.rew
            correct.rew <- correctiefactor.rew.s(p, n, alpha)
        }
        ans$cov <- ans$cov * cdelta.rew * correct.rew
        ans$call <- match.call()
        ans$method <- paste("Univariate location and scale estimation.")
        ans$alpha <- alpha
        ans$quan <- quan
        ans$raw.cov <- scale^2
        names(ans$raw.cov) <- dimn[[2]][1]
        ans$raw.center <- as.vector(center)
        names(ans$raw.center) <- dimn[[2]][1]
        ans$crit <- (1/(quan - 1)) * sum(sort((x - as.double(mcd$initmean))^2, quan)[1:quan])
        center <- ans$center
        scale <- as.vector(sqrt(ans$cov))
        ans$mcd.wt <- rep(NA, length(na.x))
        weights <- ifelse(((x - center)/scale)^2 < qchisq(0.975, p), 1, 0)

        ans$mcd.wt[ok] <- weights
        if(length(dimn[[1]]))
            names(ans$mcd.wt) <- dimn[[1]]
        ans$wt <- NULL
        if(print.it) 
        {
            cat(ans$method, "\n")
        }
        ans$X <- x
        if(length(dimn[[1]]))
            dimnames(ans$X)[[1]] <- names(ans$mcd.wt)[ok]
        else 
        {
            xx <- seq(1, length(na.x))
            dimnames(ans$X) <- list(NULL, NULL)
            dimnames(ans$X)[[1]] <- xx[ok]
        }
        class(ans) <- "mcd"
        attr(ans, "call") <- sys.call()
        return(ans)
    } #end p=1

    msg <- paste("Minimum Covariance Determinant Estimator")
    
    # If not all observations are in general position, i.e. more than h observations lie on
    # a hyperplane, the program still yields the MCD location and scatter matrix, 
    # the latter being singular (as it should be), as well as the equation of the hyperplane.
    if(mcd$exactfit != 0) 
    {
        dim(mcd$coeff) <- c(5, p)
        ans <- list()   
        ans$cov <- mcd$initcovariance
        dimnames(ans$cov) <- list(dimn[[2]], dimn[[2]])
        ans$center <- as.vector(mcd$initmean)
        if(length(dimn[[2]]))
            names(ans$center) <- dimn[[2]]
        ans$n.obs <- n
        ans$call <- match.call()
        ans$method <- msg
        if(mcd$exactfit == -1) {
            stop("The program allows for at most ", mcd$kount, 
                " observations.")
        }
        if(mcd$exactfit == -2) {
            stop("The program allows for at most ", mcd$kount, 
                " variables.")
        }
        if(mcd$exactfit == 1) {
            ans$method <- paste(ans$method, 
                "\nThe covariance matrix of the data is singular."
                )
            if(!print.it) {
                cat("The covariance matrix of the data is singular.\n"
                  )
            }
        }
        if(mcd$exactfit == 2) {
            ans$method <- paste(ans$method, 
                "\nThe covariance matrix has become singular during\nthe iterations of the MCD algorithm."
                )
            if(!print.it) {
                cat("The covariance matrix has become singular during\nthe iterations of the MCD algorithm.\n"
                  )
            }
        }
        if(p == 2) {
            ans$method <- paste(ans$method, "\nThere are", mcd$
                kount, 
                "observations in the entire dataset of\n", n, 
                "observations that lie on the line with equation\n",
                round(mcd$coeff[1, 1], digits = 4), 
                "(x_i1-m_1)+", round(mcd$coeff[1, 2], digits = 
                4), 
                "(x_i2-m_2)=0 \nwith (m_1,m_2) the mean of these observations."
                )
            if(!print.it) {
                cat("There are", mcd$kount, 
                  "observations in the entire dataset of\n", n, 
                  "observations that lie on the line with equation\n",
                  round(mcd$coeff[1, 1], digits = 4), 
                  "(x_i1-m_1)+", round(mcd$coeff[1, 2], digits
                   = 4), 
                  "(x_i2-m_2)=0 \nwith (m_1,m_2) the mean of these observations.\n"
                  )
            }
        }
        if(p == 3) {
            ans$method <- paste(ans$method, "\nThere are", mcd$
                kount, 
                "observations in the entire dataset of\n", n, 
                "observations that lie on the plane with equation \n",
                round(mcd$coeff[1, 1], digits = 4), 
                "(x_i1-m_1)+", round(mcd$coeff[1, 2], digits = 
                4), "(x_i2-m_2)+", round(mcd$coeff[1, 3], 
                digits = 4), 
                "(x_i3-m_3)=0 \nwith (m_1,m_2) the mean of these observations."
                )
            if(!print.it) {
                cat("There are", mcd$kount, 
                  "observations in the entire dataset of\n", n, 
                  "observations that lie on the plane with equation \n",
                  round(mcd$coeff[1, 1], digits = 4), 
                  "(x_i1-m_1)+", round(mcd$coeff[1, 2], digits
                   = 4), "(x_i2-m_2)+", round(mcd$coeff[1, 3], 
                  digits = 4), 
                  "(x_i3-m_3)=0 \nwith (m_1,m_2) the mean of these observations.\n"
                  )
            }
        }
        if(p > 3) {
            ans$method <- paste(ans$method, "\nThere are", mcd$
                kount, 
                " observations in the entire dataset of\n", n, 
                "observations that lie on the hyperplane with equation \na_1*(x_i1-m_1)+...+a_p*(x_ip-m_p)=0 \nwith (m_1,...,m_p) the mean\nof these observations and coefficients a_i equal to: \n"
                )
            if(!print.it) {
                cat("There are", mcd$kount, 
                  " observations in the entire dataset of\n", n,
                  "observations that lie on the hyperplane with equation \na_1*(x_i1-m_1)+...+a_p*(x_ip-m_p)=0 \nwith (m_1,...,m_p) the mean\nof these observations and coefficients a_i equal to: \n"
                  )
            }
        }
        if(p > 3) {
            for(i in 1:p) {
                ans$method <- paste(ans$method, round(mcd$coeff[
                  1, i], digits = 4))
            }
            if(!print.it)
                print(round(mcd$coeff[1,  ], digits = 4))
        }
        ans$alpha <- alpha
        ans$quan <- quan
        ans$raw.cov <- mcd$initcovariance
        dimnames(ans$raw.cov) <- list(dimn[[2]], dimn[[2]])
        ans$raw.center <- as.vector(mcd$initmean)
        if(length(dimn[[2]]))
            names(ans$raw.center) <- dimn[[2]]
        ans$crit <- 0
        ans$mcd.wt <- rep(NA, length(na.x))
        ans$mcd.wt[ok] <- mcd$weights
        if(length(dimn[[1]]))
            names(ans$mcd.wt) <- dimn[[1]]
        ans$wt <- NULL
        ans$X <- x  
        if(length(dimn[[1]]))
            dimnames(ans$X)[[1]] <- names(ans$mcd.wt)[ok]
        else {
            xx <- seq(1, length(na.x))
            dimnames(ans$X) <- list(NULL, NULL)
            dimnames(ans$X)[[1]] <- xx[ok]
        }
        
        if(print.it) {
            cat(ans$method, "\n")
        }
        class(ans) <- "mcd"
        attr(ans, "call") <- sys.call()
        return(ans)
    } #end exact fit

    weights <- mcd$weights
    weights <- as.vector(weights)   

    # Compute and apply the consistency correction factor for the reweighted cov
    if(sum(weights) == n){
        cdelta.rew <- 1
        correct.rew <- 1
    }
    else {
        qdelta.rew <- qchisq(sum(weights)/n, p)
        cdeltainvers.rew <- pgamma(qdelta.rew/2, p/2 + 1)/(sum(weights)/n)
        cdelta.rew <- 1/cdeltainvers.rew
        correct.rew <- correctiefactor.rew.s(p, n, alpha)
    }

    ans <- cov.wt(x, wt = weights, cor)
    ans$call <- match.call()
    ans$cov <- sum(weights)/(sum(weights) - 1) * ans$cov
    ans$cov <- ans$cov * cdelta.rew * correct.rew
    ans$call <- match.call()
    ans$method <- msg
    
    #vt:: add also the best found subsample to the result list
    ans$best <- sort(as.vector(mcd$best))

    # Check if the reweighted scatter matrix is singular. 
    if( - (determinant(ans$cov, log = TRUE)$modulus[1] - 0)/p > 50) {
        ans$method <- paste(ans$method, "\nThe reweighted MCD scatter matrix is singular.")
        if(!print.it) {
            cat("The reweighted MCD scatter matrix is singular.\n")
        }
        ans$alpha <- alpha
        ans$quan <- quan
        ans$raw.cov <- mcd$initcovariance
        dimnames(ans$raw.cov) <- list(dimn[[2]], dimn[[2]])
        ans$raw.center <- as.vector(mcd$initmean)
        ans$raw.weights <- weights
        ans$raw.mah <- ans$mah <- mahalanobis(x,ans$raw.center,ans$raw.cov, tol.inv = tol)
        if(length(dimn[[2]]))
            names(ans$raw.center) <- dimn[[2]]
        ans$crit <- mcd$mcdestimate    
        ans$mcd.wt <- rep(NA, length(na.x))
        ans$mcd.wt[ok] <- weights
        if(length(dimn[[1]]))
            names(ans$mcd.wt) <- dimn[[1]]
        ans$wt <- NULL
        ans$X <- x  
        if(length(dimn[[1]]))
            dimnames(ans$X)[[1]] <- names(ans$mcd.wt)[ok]
        else {
            xx <- seq(1, length(na.x))
            dimnames(ans$X) <- list(NULL, NULL)
            dimnames(ans$X)[[1]] <- xx[ok]
        }
        if(print.it)
            cat(ans$method, "\n")
        class(ans) <- "mcd"
        attr(ans, "call") <- sys.call()
        return(ans)
    }
    else {
        ans$alpha <- alpha
        ans$quan <- quan
        ans$raw.cov <- mcd$initcovariance
        dimnames(ans$raw.cov) <- list(dimn[[2]], dimn[[2]])
        ans$raw.center <- as.vector(mcd$initmean)
        ans$raw.mah <- mahalanobis(x,ans$raw.center,ans$raw.cov, tol.inv = tol)
        ans$raw.weights <- weights
        if(length(dimn[[2]]))
            names(ans$raw.center) <- dimn[[2]]
        ans$crit <- mcd$mcdestimate    
        mah <- mahalanobis(x, ans$center, ans$cov, tol.inv = tol)
        ans$mah <- mah
        weights<- ifelse(mah< qchisq(0.975, p), 1, 0)
    }
    ans$mcd.wt <- rep(NA, length(na.x))
    ans$mcd.wt[ok] <-weights
    if(length(dimn[[1]]))
        names(ans$mcd.wt) <- dimn[[1]]
    ans$wt <- NULL
    ans$X <- x
    if(length(dimn[[1]]))
        dimnames(ans$X)[[1]] <- names(ans$mcd.wt)[ok]
    else {
        xx <- seq(1, length(na.x))
        dimnames(ans$X) <- list(NULL, NULL)
        dimnames(ans$X)[[1]] <- xx[ok]
    }
    if(print.it)
        cat(ans$method, "\n")
    class(ans) <- "mcd"
    attr(ans, "call") <- sys.call()
    return(ans)
}
##  rrcov : Scalable Robust Estimators with High Breakdown Point
##
##  This program is free software; you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation; either version 2 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program; if not, write to the Free Software
##  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
##
##  I would like to thank Peter Filtzmoser for providing the initial code of 
##  some of these functions.
##

plot.mcd <- function(x,
                   which=c("all", "dd","distance","qqchi2","tolellipse"),
                   classic=FALSE,
                   ask=(which=="all" && dev.interactive()),
                   cutoff, 
                   id.n,
                   tol.inv = 1e-7, ...){
    if (!inherits(x, "mcd"))
        stop("Use only with 'mcd' objects") 
        
    covPlot(x$X, which=which, classic=classic, ask=ask, cutoff=cutoff, mcd=x, id.n=id.n, tol.inv=tol.inv, ...)
}

covPlot <- function(x, 
                   which=c("all", "dd","distance","qqchi2","tolellipse"),
                   classic=FALSE,
                   ask=FALSE,
                   mcd, 
                   cutoff, 
                   id.n,
                   tol.inv = 1e-7, ...){

##@bdescr
##  Make plots based on the covariance structure of a data set:
##  dd       -  distance-distance plot: Robust distances versus 
##              Mahalanobis distances
##  distance -  a plot of the robust distances
##  qqchi2   -  a qq-plot of the robust distances versus the 
##              quantiles of the chi-squared distribution
##  tolellipse- a tolerance ellipse
##
## Distance Plot: 
## Draw a Distance-Distance Plot: Plots the robust distances
## versus the classical Mahalanobis distances as introduced by
## Rousseeuw, P. J., and van Zomeren, B. C. (1990). Unmasking
## Multivariate Outliers and Leverage Points. Journal of the American
## Statistical Association, 85, 633-639.
##
## The dashed line is the set of points where the robust distance is 
## equal to the classical distance.
## The horizontal and vertical dotted lines are drawn at values equal cutoff
## which defaults to square root of the 97.5% quantile of a chi-squared 
## distribution with p degrees of freedom. Points beyond these lines can 
## be considered outliers.
## 
##@edescr
##
##@in  x                 : [matrix] A data.frame or matrix, n > 2*p
##@in  which          : [character] A plot option, one of:
##                            classic: index plot of the classical mahalanobis distances
##                            robust:  index plot of the robust mahalanobis distances
##                            dd:      distance-distance plot
##                            index:   parallel index plot of classical and robust distances
##                            all:     all three plots
##                          default is "all"
##@in  classic           : [logical] If true the classical plot will be displayed too
##                                   default is classic=FALSE 
##@in  mcd               : [mcd object] An object of type mcd - its attributes 
##                                      center and cov will be used
##@in  cutoff            : [number] The cutoff value for the distances 
##@in  id.n               : [number] number of observations to be identified with a label.
##                                  Defaults to the number of observations with distance
##                                  larger than cutoff 
##@in  tol.inv           : [number] tolerance to be used for computing the inverse - see 'solve'.
##                                  defaults to 1e-7 

# NOTE: The default tolerance 1e-7, will not work for some example 
#       data sets, like milk or aircraft


mydistplot <- function(x, cutoff, classic = FALSE, id.n){
##  Index Plot:
##  Plot the vector x (robust or mahalanobis distances) against 
##  the observation indexes. Identify by a label the id.n
##  observations with largest value of x. If id.n is not supplied, 
##  calculate it as the number of observations larger than cutoff.
##  Use cutoff to draw a horisontal line.
##  Use classic=FALSE/TRUE to choose the label of the vertical axes

    n <- length(x)
    if(missing(id.n))
        id.n <- length(which(x>cutoff))
    if(classic)
        ylab="Square Root of Mahalanobis distance"
    else
        ylab="Square Root of Robust distance"
    plot(x, ylab=ylab, xlab="Index", type="p")
    label(1:n, x, id.n)
    abline(h=cutoff)

    title(main="Distance Plot")
}

myddplot <- function(md, rd, cutoff, id.n){
##  Distance-Distance Plot:
##  Plot the vector y=rd (robust distances) against 
##  x=md (mahalanobis distances). Identify by a label the id.n
##  observations with largest rd. If id.n is not supplied, calculate
##  it as the number of observations larger than cutoff. Use cutoff
##  to draw a horisontal and a vertical line. Draw also a dotted line
##  with a slope 1.
    n <- length(md)
    if(missing(id.n))
        id.n <- length(which(rd>cutoff))
    xlab <- "Mahalanobis distance"
    ylab <- "Robust distance"
    plot(md, rd, xlab=xlab, ylab=ylab, type="p")
    label(md,rd,id.n)
    abline(0, 1, lty=2)
    abline(v=cutoff)
    abline(h=cutoff)

    title(main="Distance-Distance Plot")
}

qqplot <- function(x, p, cutoff, classic=FALSE, id.n){
##  Chisquare QQ-Plot:
##  Plot the vector x (robust or mahalanobis distances) against 
##  the square root of the quantiles of the chi-squared distribution
##  with p degrees of freedom.
##  Identify by a label the id.n observations with largest value of x.
##  If id.n is not supplied, calculate it as the number of observations
##  larger than cutoff.
##  Use classic=FALSE/TRUE to choose the label of the vertical axes


    ##  parameters and preconditions     

    n <- length(x)

    if(missing(cutoff))
        cutoff <- sqrt(qchisq(0.975, p))

    if(missing(id.n))
        id.n <- length(which(x>cutoff))

    qq <- sqrt(qchisq(((1:n)-1/3)/(n+1/3), p))

    x <- sort(x, index.return=TRUE)
    ix <- x$ix
    x <- x$x

    if(classic)
        ylab="Mahalanobis distance"
    else
        ylab="Robust distance"

    plot(qq, x, xlab="Square root of the quantiles of the chi-squared distribution", ylab=ylab, type="p")
    if(id.n > 0){
        ind <- (n-id.n+1):n
        xrange <- par("usr")
        xrange <- xrange[2] - xrange[1]
        text(qq[ind] + xrange/50, x[ind], ix[ind])
    }
    abline(0, 1, lty=2)
    title(main="Chisquare QQ-Plot")
}

label <- function(x, y, id.n=3){
    if(id.n > 0) {
        xrange <- par("usr")
        xrange <- xrange[2] - xrange[1]
        n <- length(y)
        ind <- sort(y, index.return=TRUE)$ix
        ind <- ind[(n-id.n+1):n]
        text(x[ind] + xrange/50, y[ind], ind)
    }
}
    
    ##  parameters and preconditions     

    if(is.vector(x) || is.matrix(x)) {
        if(!is.numeric(x))
            stop(message = "x is not a numeric dataframe or matrix.")
    }else if(is.data.frame(x)) {
        if(!all(sapply(x,data.class) == "numeric"))
            stop(message = "x is not a numeric dataframe or matrix.")
    }
    
    n <- dim(x)[1]
    p <- dim(x)[2]

    if(missing(cutoff))
        cutoff <- sqrt(qchisq(0.975, p))

    if(!missing(id.n) && !is.null(id.n)){
        id.n <- as.integer(id.n)
        if(id.n < 0 || id.n > n)
            stop("`id.n' must be in {1,..,",n,"}") 
    }

    if(missing(mcd))
        mcd <- covMcd(x)

    if(length(mcd$center)  == 0 ||  length(mcd$cov) == 0)
        stop(message = "Invalid mcd object: attributes center and cov missing!")

    if(length(mcd$center)  != p)
        stop(message = "Data set and provided center have different dimensions!")

    md <- mahalanobis(x, apply(x,2,mean), var(x), tol.inv=tol.inv)
    md <- sqrt(md)

    which <- match.arg(which)
    rd <- mahalanobis(x, mcd$center, mcd$cov, tol.inv=tol.inv)
    rd <- sqrt(rd)

    if(!classic || which == "dd")    
        par(mfrow=c(1,1), pty="m")
    else
        par(mfrow=c(1,2), pty="m")

    if (ask) {
        op <- par(ask = TRUE)
        on.exit(par(op))
    } 

    if(which == "all" || which == "distance"){    
        mydistplot(rd, cutoff, id.n=id.n)                     # index plot of mahalanobis distances
        if(classic)
            mydistplot(md, cutoff, classic=TRUE, id.n=id.n)   # index plot of robust distances
    }

    if(which == "all" || which == "dd"){    
        myddplot(md, rd, cutoff=cutoff, id.n=id.n)    # distance-distance plot
    }

    if(which == "all" || which == "qqchi2"){    
        qqplot(rd, p, cutoff=cutoff, id.n=id.n)     # qq-plot of the robust distances versus the 
                                                    # quantiles of the chi-squared distribution
        if(classic)
            qqplot(md, p, cutoff=cutoff, classic=TRUE, id.n=id.n)
                                                    # qq-plot of the mahalanobis distances
    }

    if(which == "all" || which == "tolellipse"){    
       tolellipse(x, mcd=mcd, cutoff=cutoff, id.n=id.n, classic=classic, tol.inv=tol.inv)
                                                # qq-plot of the robust distances versus the 
                                                # quantiles of the chi-squared distribution
    }
}

ddplot <- function(x,...){
    covPlot(x, which="dd", ...)
}

distplot <- function(x,...){
    covPlot(x, which="distance", ...)
}

chi2qqplot <- function(x,...){
    covPlot(x, which="qqchi2", ...)
}

ellipse <- function(x,...){
    covPlot(x, which="tolellipse", ...)
}
##  rrcov : Scalable Robust Estimators with High Breakdown Point
##
##  This program is free software; you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation; either version 2 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program; if not, write to the Free Software
##  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
##

plot.lts <- function(x,
                   which=c("all", "rqq","rindex", "rfit", "rdiag"),
                   classic=FALSE,
                   ask=(which=="all" && dev.interactive()),
                   id.n=3, ...){
    if (!inherits(x, "lts"))
        stop("Use only with 'lts' objects") 
        
    ltsPlot(x, which, classic, ask, id.n, ...)
}

ltsPlot <- function(x, 
                   which = c("all", "rqq","rindex", "rfit", "rdiag"),
                   classic = FALSE,
                   ask=FALSE,
                   id.n=3, ...){
##@bdescr
##  Make plots for model checking and outlier detection based on 
##      the LTS regression estimates:
##  rqq      -  normal quantile plot of the LTS and LS residuals
##  rindex   -  standardized LTS/LS Residuals versus index
##  rfit     -  standardized LTS/LS Residuals versus fitted values
##  rdiag    -  regression diagnostic plot
##
##@edescr
##
##@in  x               : [object] An lts object
##@in  which          : [character] A plot option, one of:
##                            rqq:    
##                            rdiag:  
##                            rfit:   
##                            rindex: 
##                          default is "rqq"
##@in  classic           : [logical] If true the classical plot will be displayed too
##                                   default is classic=FALSE 
##@in  id.n               : [number] number of observations to be identified with a label.
##                                  Defaults to 3

label <- function(x, y, ord, lab, id.n, ...)
{
    if(id.n) {
        n <- length(y)
        which <- order(ord)[(n - id.n + 1):n]
        if(missing(lab))
            lab <- which
        else
            lab <- lab[which]
        # how to adjust the labels?
        # a) adj=0.1
        # b) x=x+xrange
        # c) pos=4 (to the left of the observation)
        # d) additionaly to pos specify offset=0.2 (fraction of a character)
        xrange <- par("usr")
        xrange <- (xrange[2] - xrange[1])/50
        text(x[which], y[which], pos=4, offset=0.2, lab, ...)
    }
}

# The R function 'qqline' (package::stats) adds a line to a 
# normal quantile-quantile plot which passes through the 
# first and third quartiles. In S this function returns the
# slope and intercept of the line, but not in R. 
# Here we need the slope and intercept in order to sort the
# residuals according to their distance from the line.

myqqline <- function(y, datax = FALSE, ...){
    y <- quantile(y[!is.na(y)],c(0.25, 0.75))
    x <- qnorm(c(0.25, 0.75))
    if(datax) {
        slope <- diff(x)/diff(y)
        int <- x[1] - slope*y[1]
    } else {
        slope <- diff(y)/diff(x)
        int <- y[1]-slope*x[1]
    }
    abline(int, slope, ...)
    invisible(list(int=int, slope=slope))
} 

myqqplot <- function(r, classic = FALSE, lab, id.n, ...){
##  Normal QQ-plot of residuals:
##  Produces a Quantile-Quantile plot in which the vector r is plotted 
##  against the quantiles of a standard normal distribution.
##

    mgp = c(2.5, 1, 0)  # set the margin line (in 'mex' units) for the:
                        # - axis title, 
                        # - axis labels and 
                        # - axis line. 
                        # The default is 'c(3, 1, 0)'.
    xlab <- "Quantiles of the standard normal distribution"
    ylab <- "Standardized LTS residual"
    if(classic)
        ylab <- "Standardized LS residual"

    qq <- qqnorm(r, mgp = mgp, xlab=xlab, ylab=ylab, ...)
    ll <- myqqline(r, lty = 2, ...)
    ord <- abs(qq$y - ll$int - ll$slope * qq$x)
    label(qq$x, qq$y, ord, lab, id.n, ...)
    title(main="Normal Q-Q plot")
}

indexplot <- function(r, scale, classic = FALSE, lab, id.n, ...){
##  Index plot of standardized residuals:
##  Plot the vector r (LTS or LS residuals) against 
##  the observation indexes. Identify by a label the id.n
##  observations with largest value of r. 
##  Use classic=FALSE/TRUE to choose the label of the vertical axes

    mgp = c(2.5, 1, 0)  # set the margin line (in 'mex' units) for the:
                        # - axis title, 
                        # - axis labels and 
                        # - axis line. 
                        # The default is 'c(3, 1, 0)'.
    xlab <- "Index"
    ylab <- "Standardized LTS residual"
    if(classic)
        ylab <- "Standardized LS residual"

    x <- 1:length(r)
    y <- r/scale
    ylim <- c(min(-3, min(y)), max(3, max(y)))

    plot(x, y, ylim=ylim, mgp = mgp, xlab=xlab, ylab=ylab, ...)
    label(x, y, ord=abs(y), lab, id.n, ...)
    abline(h = -2.5, ...)
    abline(h = 0, lty = 4, ...)
    abline(h = 2.5, ...)

    mtext("-2.5", side = 4, line = 1.2, at = -2.5, ...)
    mtext("2.5", side = 4, line = 1.2, at = 2.5, ...)

    title(main="Residuals vs Index")
}

fitplot <- function(obj, classic = FALSE, lab, id.n, ...){
##  Standardized residuals vs Fitted values plot:
##  Plot the vector r (LTS or LS residuals) against 
##  the corresponding fitted values. Identify by a 
##  label the id.n observations with largest value of r. 
##  Use classic=FALSE/TRUE to choose the label of the vertical axes

    mgp = c(2.5, 1, 0)  # set the margin line (in 'mex' units) for the:
                        # - axis title, 
                        # - axis labels and 
                        # - axis line. 
                        # The default is 'c(3, 1, 0)'.

#    x <- obj$X %*% as.matrix(obj$coef)
    x <- obj$fitted.values
    y <- obj$residuals/obj$scale
    ylim <- c(min(-3, min(y)), max(3, max(y)))
    yname <- names(obj$scale)
    xlab <- paste("Fitted :", yname)
    ylab <- "Standardized LTS residual"
    if(classic)
        ylab <- "Standardized LS residual"

    plot(x, y, ylim=ylim, mgp = mgp, xlab=xlab, ylab=ylab, ...)
    label(x, y, ord=abs(y), lab, id.n, ...)
    abline(h = -2.5, ...)
    abline(h = 0, lty = 4, ...)
    abline(h = 2.5, ...)

    mtext("-2.5", side = 4, line = 1.2, at = -2.5, ...)
    mtext("2.5", side = 4, line = 1.2, at = 2.5, ...)
    
    title(main="Residuals vs Fitted")
    
}


rdiag <- function(obj, classic = FALSE, lab, id.n, ...){
##  Standardized residuals vs Fitted values plot:
##  Plot the vector r (LTS or LS residuals) against 
##  the corresponding fitted values. Identify by a 
##  label the id.n observations with largest value of r. 
##  Use classic=FALSE/TRUE to choose the label of the vertical axes

    p <- if(obj$intercept) length(obj$coef) - 1 else length(obj$coef)
    if(p <= 0)
        warning("Diagnostic plot is not available for univar\niate location and scale estimation")

    if(is.null(obj$RD))
        stop("option mcd=F was set in ltsreg.")
    if(obj$RD[1] == "singularity")
        stop("The MCD covariance matrix was singular.")

    mgp = c(2.5, 1, 0)  # set the margin line (in 'mex' units) for the:
                        # - axis title, 
                        # - axis labels and 
                        # - axis line. 
                        # The default is 'c(3, 1, 0)'.

    xlab <- "Robust distance computed by MCD"
    ylab <- "Standardized LTS residual"
    if(classic){
        xlab <- "Mahalanobis distance"
        ylab <- "Standardized LS residual"
    }

    quant <- max(c(sqrt(qchisq(0.975, p)), 2.5))
    x <- obj$RD
    y <- obj$residuals/obj$scale
    xlim <- c(0, max(quant + 0.1, max(x)))
    ylim <- c(min(-3, min(y)), max(3, max(y)))

    plot(x, y, ylim=ylim, mgp = mgp, xlab=xlab, ylab=ylab, ...)
    ord = apply(abs(cbind(x/2.5, y/quant)), 1, max)
    label(x, y, ord=ord, lab, id.n, ...)
    abline(h = -2.5, ...)
    abline(h = 2.5, ...)
    abline(v = quant, ...)

    mtext("-2.5", side = 4, line = 1.2, at = -2.5, ...)
    mtext("2.5", side = 4, line = 1.2, at = 2.5,...)

    title(main="Regression Diagnostic Plot")
}

    ##  parameters and preconditions     

    which <- match.arg(which)
    r <- residuals(x)
    n <- length(r) 
    if(!missing(id.n) && !is.null(id.n)){
        id.n <- as.integer(id.n)
        if(id.n < 0 || id.n > n)
            stop("`id.n' must be in {1,..,",n,"}") 
    }

    if(!classic)
        par(mfrow=c(1,1), pty="m")
    else {
        par(mfrow=c(1,2), pty="m")
        
        # calculate the LS regression (using LTS with alpha = 1)
        # if intercept, obj$X is augmented with a column of 1s - remove it

        if(x$intercept &&                         # model with intercept
           length(dim(x$X)) == 2 &&               # X is 2-dimensional
           dim(x$X)[2] > 1 &&                     # X has more than 1 column
           all(x$X[,dim(x$X)[2]]==1))             # the last column of X is all 1s
            X <- x$X[,1:(dim(x$X)[2]-1)]
        else
            X <- x$X
        obj.cl <- ltsReg(X, x$Y, intercept=x$intercept, alpha=1)
    }

    if (ask) {
        op <- par(ask = TRUE)
        on.exit(par(op))
    } 
    
    if(which == "all" || which == "rqq"){    
        myqqplot(x$residuals, id.n=id.n, ...)                                   # normal QQ-plot of the LTS residuals
        if(classic)
            myqqplot(obj.cl$residuals, classic=TRUE, id.n=id.n, ...)            # normal QQ-plot of the LS residuals
    }

    if(which == "all" || which == "rindex"){    
        indexplot(x$residuals, x$scale, id.n=id.n, ...)                       # index plot of the LTS residuals
        if(classic)
            indexplot(obj.cl$residuals, obj.cl$scale, classic=TRUE, id.n=id.n, ...)     # index plot of the LS residuals
    }

    if(which == "all" || which == "rfit"){    
        fitplot(x, id.n=id.n, ...)                       
        if(classic)
            fitplot(obj.cl, classic=TRUE, id.n=id.n, ...)     
    }

    if(which == "all" || which == "rdiag"){    
        rdiag(x, id.n=id.n, ...)                       
        if(classic)
            rdiag(obj.cl, classic=TRUE, id.n=id.n, ...)     
    }
}
##  rrcov : Scalable Robust Estimators with High Breakdown Point
##
##  This program is free software; you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation; either version 2 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program; if not, write to the Free Software
##  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
##
##  I would like to thank Peter Rousseeuw and Katrien van Driessen for 
##  providing the initial code of this function.


ltsReg <- function (x, y, 
                    intercept=TRUE, 
                    alpha=NULL, 
                    nsamp=500, 
                    adjust=FALSE, 
                    mcd=TRUE, 
                    qr.out=FALSE, 
                    yname=NULL, 
                    seed=0) 
{

    quan.f <- function(alpha, n, rk) {
        quan <- floor(2*floor((n+rk+1)/2) - n + 2*(n - floor((n+rk+1)/2)) * alpha)
        return(quan)
    }

    correctiefactor.s <- function(p, intercept = intercept, n, alpha) {
        if (intercept == TRUE) {
            p <- p - 1
        }
        if (p == 0) {
            fp.500.n <- 1 - exp(0.262024211897096) * 1/n^{
                0.604756680630497
            }
            fp.875.n <- 1 - exp(-0.351584646688712) * 1/n^{
                1.01646567502486
            }
            if ((0.5 <= alpha) && (alpha <= 0.875)) {
                fp.alpha.n <- fp.500.n + (fp.875.n - fp.500.n)/0.375 * 
                  (alpha - 0.5)
                fp.alpha.n <- sqrt(fp.alpha.n)
            }
            if ((0.875 < alpha) && (alpha < 1)) {
                fp.alpha.n <- fp.875.n + (1 - fp.875.n)/0.125 * 
                  (alpha - 0.875)
                fp.alpha.n <- sqrt(fp.alpha.n)
            }
        }
        else {
            if (p == 1) {
                if (intercept == TRUE) {
                  fp.500.n <- 1 - exp(0.630869217886906) * 1/n^{
                    0.650789250442946
                  }
                  fp.875.n <- 1 - exp(0.565065391014791) * 1/n^{
                    1.03044199012509
                  }
                }
                else {
                  fp.500.n <- 1 - exp(-0.0181777452315321) * 
                    1/n^{
                    0.697629772271099
                  }
                  fp.875.n <- 1 - exp(-0.310122738776431) * 1/n^{
                    1.06241615923172
                  }
                }
            }
            if (p > 1) {
                if (intercept == TRUE) {
                  coefgqpkwad875 <- matrix(c(-0.458580153984614, 
                    1.12236071104403, 3, -0.267178168108996, 
                    1.1022478781154, 5), ncol = 2, byrow = FALSE)
                  dimnames(coefgqpkwad875) <- list(c("alfaq", 
                    "betaq", "qwaarden"), c("coefgqpkwad875mint.q3", 
                    "coefgqpkwad875mint.q5"))
                  coefeqpkwad500 <- matrix(c(-0.746945886714663, 
                    0.56264937192689, 3, -0.535478048924724, 
                    0.543323462033445, 5), ncol = 2, byrow = FALSE)
                  dimnames(coefeqpkwad500) <- list(c("alfaq", 
                    "betaq", "qwaarden"), c("coefeqpkwad500mint.q3", 
                    "coefeqpkwad500mint.q5"))
                }
                else {
                  coefgqpkwad875 <- matrix(c(-0.251778730491252, 
                    0.883966931611758, 3, -0.146660023184295, 
                    0.86292940340761, 5), ncol = 2, byrow = FALSE)
                  dimnames(coefgqpkwad875) <- list(c("alfaq", 
                    "betaq", "qwaarden"), c("coefgqpkwad875zint.q3", 
                    "coefgqpkwad875zint.q5"))
                  coefeqpkwad500 <- matrix(c(-0.487338281979106, 
                    0.405511279418594, 3, -0.340762058011, 0.37972360544988, 
                    5), ncol = 2, byrow = FALSE)
                  dimnames(coefeqpkwad500) <- list(c("alfaq", 
                    "betaq", "qwaarden"), c("coefeqpkwad500zint.q3", 
                    "coefeqpkwad500zint.q5"))
                }
                y1.500 <- 1 + coefeqpkwad500[1, 1] * 1/p^{
                  coefeqpkwad500[2, 1]
                }
                y2.500 <- 1 + coefeqpkwad500[1, 2] * 1/p^{
                  coefeqpkwad500[2, 2]
                }
                y1.875 <- 1 + coefgqpkwad875[1, 1] * 1/p^{
                  coefgqpkwad875[2, 1]
                }
                y2.875 <- 1 + coefgqpkwad875[1, 2] * 1/p^{
                  coefgqpkwad875[2, 2]
                }
                y1.500 <- log(1 - y1.500)
                y2.500 <- log(1 - y2.500)
                y.500 <- c(y1.500, y2.500)
                A.500 <- matrix(c(1, log(1/(coefeqpkwad500[3, 
                  1] * p^2)), 1, log(1/(coefeqpkwad500[3, 2] * 
                  p^2))), ncol = 2, byrow = TRUE)
                coeffic.500 <- solve(A.500, y.500)
                y1.875 <- log(1 - y1.875)
                y2.875 <- log(1 - y2.875)
                y.875 <- c(y1.875, y2.875)
                A.875 <- matrix(c(1, log(1/(coefgqpkwad875[3, 
                  1] * p^2)), 1, log(1/(coefgqpkwad875[3, 2] * 
                  p^2))), ncol = 2, byrow = TRUE)
                coeffic.875 <- solve(A.875, y.875)
                fp.500.n <- 1 - exp(coeffic.500[1]) * 1/n^{
                  coeffic.500[2]
                }
                fp.875.n <- 1 - exp(coeffic.875[1]) * 1/n^{
                  coeffic.875[2]
                }
            }
            if ((0.5 <= alpha) && (alpha <= 0.875)) {
                fp.alpha.n <- fp.500.n + (fp.875.n - fp.500.n)/0.375 * 
                  (alpha - 0.5)
            }
            if ((0.875 < alpha) && (alpha <= 1)) {
                fp.alpha.n <- fp.875.n + (1 - fp.875.n)/0.125 * 
                  (alpha - 0.875)
            }
        }
        return(1/fp.alpha.n)
    }
    correctiefactor.rew.s <- function(p, intercept = intercept, n, alpha) {
        if (intercept == TRUE) {
            p <- p - 1
        }
        if (p == 0) {
            fp.500.n <- 1 - (exp(1.11098143415027) * 1)/n^{
                1.5182890270453
            }
            fp.875.n <- 1 - (exp(-0.66046776772861) * 1)/n^{
                0.88939595831888
            }
            if ((0.5 <= alpha) && (alpha <= 0.875)) {
                fp.alpha.n <- fp.500.n + (fp.875.n - fp.500.n)/0.375 * 
                  (alpha - 0.5)
                fp.alpha.n <- sqrt(fp.alpha.n)
            }
            if ((0.875 < alpha) && (alpha < 1)) {
                fp.alpha.n <- fp.875.n + (1 - fp.875.n)/0.125 * 
                  (alpha - 0.875)
                fp.alpha.n <- sqrt(fp.alpha.n)
            }
        }
        else {
            if (p == 1) {
                if (intercept == TRUE) {
                  fp.500.n <- 1 - (exp(1.58609654199605) * 1)/n^{
                    1.46340162526468
                  }
                  fp.875.n <- 1 - (exp(0.391653958727332) * 1)/n^{
                    1.03167487483316
                  }
                }
                else {
                  fp.500.n <- 1 - (exp(0.6329852387657) * 1)/n^{
                    1.40361879788014
                  }
                  fp.875.n <- 1 - (exp(-0.642240988645469) * 
                    1)/n^{
                    0.926325452943084
                  }
                }
            }
            if (p > 1) {
                if (intercept == TRUE) {
                  coefqpkwad875 <- matrix(c(-0.474174840843602, 
                    1.39681715704956, 3, -0.276640353112907, 
                    1.42543242287677, 5), ncol = 2, byrow = FALSE)
                  dimnames(coefqpkwad875) <- list(c("alfaq", 
                    "betaq", "qwaarden"), c("coefqpkwad875mint.q3", 
                    "coefqpkwad875mint.q5"))
                  coefqpkwad500 <- matrix(c(-0.773365715932083, 
                    2.02013996406346, 3, -0.337571678986723, 
                    2.02037467454833, 5), ncol = 2, byrow = FALSE)
                  dimnames(coefqpkwad500) <- list(c("alfaq", 
                    "betaq", "qwaarden"), c("coefqpkwad500mint.q3", 
                    "coefqpkwad500mint.q5"))
                }
                else {
                  coefqpkwad875 <- matrix(c(-0.267522855927958, 
                    1.17559984533974, 3, -0.161200683014406, 
                    1.21675019853961, 5), ncol = 2, byrow = FALSE)
                  dimnames(coefqpkwad875) <- list(c("alfaq", 
                    "betaq", "qwaarden"), c("coefqpkwad875zint.q3", 
                    "coefqpkwad875zint.q5"))
                  coefqpkwad500 <- matrix(c(-0.417574780492848, 
                    1.83958876341367, 3, -0.175753709374146, 
                    1.8313809497999, 5), ncol = 2, byrow = FALSE)
                  dimnames(coefqpkwad500) <- list(c("alfaq", 
                    "betaq", "qwaarden"), c("coefqpkwad500zint.q3", 
                    "coefqpkwad500zint.q5"))
                }
                y1.500 <- 1 + (coefqpkwad500[1, 1] * 1)/p^{
                  coefqpkwad500[2, 1]
                }
                y2.500 <- 1 + (coefqpkwad500[1, 2] * 1)/p^{
                  coefqpkwad500[2, 2]
                }
                y1.875 <- 1 + (coefqpkwad875[1, 1] * 1)/p^{
                  coefqpkwad875[2, 1]
                }
                y2.875 <- 1 + (coefqpkwad875[1, 2] * 1)/p^{
                  coefqpkwad875[2, 2]
                }
                y1.500 <- log(1 - y1.500)
                y2.500 <- log(1 - y2.500)
                y.500 <- c(y1.500, y2.500)
                A.500 <- matrix(c(1, log(1/(coefqpkwad500[3, 
                  1] * p^2)), 1, log(1/(coefqpkwad500[3, 2] * 
                  p^2))), ncol = 2, byrow = TRUE)
                coeffic.500 <- solve(A.500, y.500)
                y1.875 <- log(1 - y1.875)
                y2.875 <- log(1 - y2.875)
                y.875 <- c(y1.875, y2.875)
                A.875 <- matrix(c(1, log(1/(coefqpkwad875[3, 
                  1] * p^2)), 1, log(1/(coefqpkwad875[3, 2] * 
                  p^2))), ncol = 2, byrow = TRUE)
                coeffic.875 <- solve(A.875, y.875)
                fp.500.n <- 1 - (exp(coeffic.500[1]) * 1)/n^{
                  coeffic.500[2]
                }
                fp.875.n <- 1 - (exp(coeffic.875[1]) * 1)/n^{
                  coeffic.875[2]
                }
            }
            if ((0.5 <= alpha) && (alpha <= 0.875)) {
                fp.alpha.n <- fp.500.n + (fp.875.n - fp.500.n)/0.375 * 
                  (alpha - 0.5)
            }
            if ((0.875 < alpha) && (alpha <= 1)) {
                fp.alpha.n <- fp.875.n + (1 - fp.875.n)/0.125 * 
                  (alpha - 0.875)
            }
        }
        return(1/fp.alpha.n)
    }
    
    
    if (is.vector(y) || (is.matrix(y) && !is.data.frame(y))) {
        if (!is.numeric(y)) 
            stop(message = "y is not a numeric dataframe or vector.")
    }
    if ((!is.matrix(y) && !is.vector(y)) || is.data.frame(y)) {
        if ((!is.data.frame(y) && !is.numeric(y)) || (!all(sapply(y, data.class) == "numeric"))) 
            stop(message = "y is not a numeric dataframe or vector.")
    }
    
    y <- as.matrix(y)
    if (dim(y)[2] != 1) 
        stop(message = "y is not onedimensional.")
    
    if (missing(x)) {
        x <- rep(1, nrow(y))
        if (is.vector(x) || (is.matrix(x) && !is.data.frame(x))) {
            if (!is.numeric(x)) 
                stop(message = "x is not a numeric dataframe or matrix.")
        }
        if ((!is.matrix(x) && !is.vector(x)) || is.data.frame(x)) {
            if ((!is.data.frame(x) && !is.numeric(x)) || (!all(sapply(x, 
                data.class) == "numeric"))) 
                stop(message = "x is not a numeric dataframe or matrix.")
        }
        if (!is.matrix(x)) 
            x <- array(x, c(length(x), 1), list(names(x), deparse(substitute(x))))
        x <- as.matrix(x)
        dimny <- dimnames(y)
        dimnx <- dimnames(x)
        na.x <- !is.finite(x %*% rep(1, ncol(x)))
        na.y <- !is.finite(y)
        if (nrow(na.x) != nrow(na.y)) 
            stop("Number of observations in x and y not equal")
        ok <- !(na.x | na.y)
        y <- y[ok, , drop = FALSE]
        dy <- nrow(y)
        rownames <- dimny[[1]]
        yn <- if (!is.null(yname)) yname else dimny[[2]]
        if (!length(yn)) 
            yn <- "Y"
        storage.mode(y) <- "double"
        x <- x[ok, , drop = FALSE]
        storage.mode(x) <- "double"
        dx <- dim(x)
        if (!length(dx)) 
            stop("All observations have missing values!")
        n <- dx[1]
    }else {
        if (is.vector(x) || (is.matrix(x) && !is.data.frame(x))) {
            if (!is.numeric(x)) 
                stop(message = "x is not a numeric dataframe or matrix.")
        }
        if ((!is.matrix(x) && !is.vector(x)) || is.data.frame(x)) {
            if ((!is.data.frame(x) && !is.numeric(x)) || (!all(sapply(x, data.class) == "numeric"))) 
                stop(message = "x is not a numeric dataframe or matrix.")
        }

        #vt:: if the data is supplied as a data.frame, the following expressions results in an error
        # as workaround convert the data.frame to a matrix
        if(is.data.frame(x))
            x <- as.matrix(x)
    
        if (!is.matrix(x)) 
            x <- array(x, c(length(x), 1), list(names(x), deparse(substitute(x))))
        x <- as.matrix(x)
        dimny <- dimnames(y)
        dimnx <- dimnames(x)
        na.x <- !is.finite(x %*% rep(1, ncol(x)))
        na.y <- !is.finite(y)
        if (nrow(na.x) != nrow(na.y)) 
            stop("Number of observations in x and y not equal")
        ok <- !(na.x | na.y)
        y <- y[ok, , drop = FALSE]
        dy <- nrow(y)
        rownames <- dimny[[1]]
        yn <- if (!is.null(yname)) 
            yname
        else dimny[[2]]
        if (!length(yn)) 
            yn <- "Y"
        storage.mode(y) <- "double"
        x <- x[ok, , drop = FALSE]
        storage.mode(x) <- "double"
        dx <- dim(x)
        if (!length(dx)) 
            stop("All observations have missing values!")
        n <- dx[1]
        constantcolom <- function(x) {
            c1 <- range(x)
            c1[1] == c1[2]
        }
        if (sum(apply(x, 2, constantcolom)) > 0) 
            stop("There is at least one constant column. Remove this column and set intercept=T")
    }
    dn <- dimnames(x)
    xn <- dn[[2]]
    if (!length(xn)) 
        if (dx[2] > 1) 
            xn <- paste("X", 1:dx[2], sep = "")
        else xn <- "X"
    X <- x
    dimnames(X) <- list(NULL, xn)
    y <- as.vector(y)
    if (all(x == 1)) {
        if (length(alpha)) {
            if (alpha < 1/2) 
                stop("alpha is out of range!")
            else if (alpha > 1) 
                stop("alpha is greater than 1")
            quan <- quan.f(alpha, n, dx[2])
        }
        else {
            alpha <- 1/2
            quan <- quan.f(alpha, n, dx[2])
        }
        storage.mode(quan) <- "integer"
        initcov <- matrix(0, nrow = 1, ncol = 1)
        initmean <- matrix(0, nrow = 1, ncol = 1)
        inbest <- matrix(0, nrow = quan, ncol = 1)
        weights <- matrix(0, nrow = n, ncol = 1)
        coeff <- matrix(0, nrow = 5, ncol = 1)
        adjustcov <- matrix(0, nrow = 1, ncol = 1)
        storage.mode(initcov) <- "double"
        storage.mode(initmean) <- "double"
        storage.mode(inbest) <- "integer"
        storage.mode(weights) <- "integer"
        storage.mode(coeff) <- "double"
        storage.mode(adjustcov) <- "double"
        storage.mode(y) <- "double"
        storage.mode(seed) <- "integer"
        exactfit <- 0
        p <- 1
        xbest <- NULL
        if (alpha == 1) {
            scale <- sqrt(cov.wt(x)$cov)
            center <- as.vector(mean(x))
        }else {
            sh <- .Fortran("rffastmcd", 
                    as.matrix(y), 
                    as.integer(n), 
                    as.integer(p), 
                    as.integer(quan), 
                    nsamp = 0, 
                    initcovariance = initcov, 
                    initmean = initmean, 
                    inbest = inbest, 
                    mcdestimate = 0, 
                    weights = weights, 
                    as.integer(exactfit), 
                    coeff = coeff, 
                    kount = 0, 
                    adjustcov = adjustcov,
                    seed,
                    PACKAGE="rrcov")

            y <- as.vector(y)
            center <- as.double(sh$initmean)
            qalpha <- qchisq(quan/n, 1)
            calphainvers <- pgamma(qalpha/2, 1/2 + 1)/(quan/n)
            calpha <- 1/calphainvers
            correct <- correctiefactor.s(1, intercept = intercept, 
                n, alpha)
            scale <- sqrt(as.double(sh$initcovariance)) * sqrt(calpha) * 
                correct
            xbest <- sort(as.vector(sh$inbest))
        }
        resid <- y - center
        ans <- list()
        ans$best <- xbest
        ans$coefficients <- center
        ans$alpha <- alpha
        ans$quan <- quan
        ans$raw.resid <- resid/scale
        weights <- rep(NA, n)
        if (abs(scale) < 1e-07) {
            weights <- ifelse(abs(resid) < 1e-07, 1, 0)
            ans$scale <- ans$raw.scale <- 0
            ans$crit <- 0
            ans$coefficients <- ans$raw.coefficients <- center
        }
        if (abs(scale) >= 1e-07) {
            ans$raw.scale <- scale
            ans$raw.coefficients <- center
            quantiel <- qnorm(0.9875)
            weights <- ifelse(abs(resid/scale) <= quantiel, 1, 
                0)
            reweighting <- cov.wt(y, wt = weights)
            ans$coefficients <- reweighting$center
            ans$scale <- sqrt(sum(weights)/(sum(weights) - 1) * 
                reweighting$cov)
            resid <- y - ans$coefficients
            ans$crit <- sum(sort((y - center)^2, quan)[1:quan])
            if (sum(weights) == n) {
                cdelta.rew <- 1
                correct.rew <- 1
            }
            else {
                qdelta.rew <- qchisq(sum(weights)/n, 1)
                cdeltainvers.rew <- pgamma(qdelta.rew/2, 1/2 + 
                  1)/(sum(weights)/n)
                cdelta.rew <- sqrt(1/cdeltainvers.rew)
                correct.rew <- correctiefactor.rew.s(1, intercept = intercept, 
                  n, alpha)
            }
            ans$scale <- ans$scale * cdelta.rew * correct.rew
            quantiel <- qnorm(0.9875)
            weights <- ifelse(abs(resid/ans$scale) <= quantiel, 
                1, 0)
        }
        ans$resid <- resid/ans$scale
        ans$rsquared <- 0
        ans$residuals <- rep(NA, length(na.y))
        ans$residuals[ok] <- resid
        ans$lts.wt <- rep(NA, length(na.y))
        ans$lts.wt[ok] <- weights
        ans$intercept <- intercept
        ans$method <- paste("Univariate location and scale estimation.")
        if (abs(scale) < 1e-07) 
            ans$method <- paste(ans$method, "\nMore than half of the data are equal!")
        names(ans$coefficients) <- names(ans$raw.coefficients) <- yn
        names(ans$scale) <- names(ans$raw.scale) <- yn
        names(ans$rsquared) <- yn
        names(ans$crit) <- yn
        names(ans$residuals) <- rownames
        names(ans$lts.wt) <- rownames
        ans$X <- x
        ans$Y <- y                  # VT:: 01.09.2004 - add y to the result object
        if (length(rownames)) 
            dimnames(ans$X)[[1]] <- rownames[ok]
        else {
            xx <- seq(1, length(na.x))
            dimnames(ans$X) <- list(NULL, NULL)
            dimnames(ans$X)[[1]] <- xx[ok]
        }
        class(ans) <- "lts"
        attr(ans, "call") <- sys.call()
        return(ans)
    }
    ans <- list()
    if (intercept) {
        dx <- dx + c(0, 1)
        xn <- c(xn, "Intercept")
        x <- array(c(x, rep(1, n)), dx, dimnames = list(dn[[1]], 
            xn))
    }
    p <- dx[2]
    if (n <= 2 * p) 
        stop("Need more than twice as many observations as variables.")
    if (length(alpha)) {
        if (alpha > 1) 
            stop("alpha is greater than 1")
        if (alpha == 1) {                                   # alpha == 1 -----------------------
            z <- lsfit(x, y, intercept = FALSE)
            ans$raw.coefficients[2:p] <- z$coef[1:(p - 1)]
            ans$raw.coefficients[1] <- z$coef[p]
            ans$alpha <- alpha
            ans$quan <- quan <- n           # VT:: 01.09.2004 - bug in alpha=1 
                                            # (ans$quan was not set)
            names(ans$raw.coefficients)[2:p] <- xn[1:(p - 1)]
            names(ans$raw.coefficients)[1] <- xn[p]
            s0 <- sqrt((1/(n - p)) * sum(z$residuals^2))
            weights <- rep(NA, n)
            if(abs(s0) < 1e-07) {
                fitted <- x %*% z$coef
                weights <- ifelse(abs(z$residuals) <= 1e-07, 1, 0)
                ans$scale <- ans$raw.scale <- 0
                ans$coefficients <- ans$raw.coefficients
            }
            else {
                ans$raw.scale <- s0
                ans$raw.resid <- ans$residuals/ans$raw.scale
                weights <- ifelse(abs(z$residuals/s0) <= qnorm(0.9875), 1, 0)
                
                # vt:: weights has to be a vector instead of a matrix - 
                #      to avoid "Error in x * wtmult : non-conformable arrays"
                # 
                weights <- as.vector(weights)
                z <- lsfit(x, y, wt = weights, intercept = FALSE)
                ans$coefficients[2:p] <- z$coef[1:(p - 1)]
                ans$coefficients[1] <- z$coef[p]
                fitted <- x %*% z$coef
                ans$scale <- sqrt(sum(weights * z$residuals^2)/(sum(weights) - 
                  1))
                if (sum(weights) == n) {
                  cdelta.rew <- 1
                }
                else {
                  cdelta.rew <- (1/sqrt(1 - ((2 * n)/(sum(weights) * 
                    (1/qnorm((sum(weights) + n)/(2 * n))))) * 
                    dnorm(1/(1/(qnorm((sum(weights) + n)/(2 * 
                      n)))))))
                }
                ans$scale <- ans$scale * cdelta.rew
                weights <- ifelse(abs(z$residuals/ans$scale) <= 
                  qnorm(0.9875), 1, 0)
                ans$resid <- z$residuals/ans$scale
            }
            names(ans$coefficients)[2:p] <- xn[1:(p - 1)]
            names(ans$coefficients)[1] <- xn[p]
            ans$crit <- sum(z$residuals^2)
            if (intercept) {
                s1 <- sum(z$residuals^2)
                center <- mean(y)
                sh <- sum((y - center)^2)
                ans$rsquared <- 1 - (s1/sh)
            }
            else {
                s1 <- sum(z$residuals^2)
                sh <- sum(y^2)
                ans$rsquared <- 1 - (s1/sh)
            }
            if (ans$rsquared > 1) {
                ans$rsquared <- 1
            }
            if (ans$rsquared < 0) {
                ans$rsquared <- 0
            }
            ans$residuals <- rep(NA, length(na.y))
            ans$residuals[ok] <- z$residuals
            ans$lts.wt <- matrix(NA, length(na.y))
            ans$lts.wt[ok] <- weights
            ans$intercept <- intercept
            ans$method <- paste("Least Squares Regression.")
            if (abs(s0) < 1e-07) 
                ans$method <- paste(ans$method, , "\nAn exact fit was found!")
            if (mcd) {
                # vt:: changed name of the function
                # mcd <- cov.mcd.default(X, print.it = FALSE, alpha = 1)
                mcd <- covMcd(X, print.it = FALSE, alpha = 1)
                if(-(determinant(mcd$cov, log = TRUE)$modulus[1])/p > 50) {
                  ans$RD[1] <- "singularity"
                }else {
                  ans$RD <- rep(NA, length(na.y))
                  ans$RD[ok] <- sqrt(mahalanobis(X, mcd$center, mcd$cov))
                  names(ans$RD) <- rownames
                }
            }
            names(ans$residuals) <- rownames
            names(ans$lts.wt) <- rownames
            names(ans$scale) <- names(ans$raw.scale) <- yn
            names(ans$rsquared) <- yn
            names(ans$crit) <- yn
            ans$X <- x
            ans$Y <- y          # VT:: 01.09.2004 - add y to the result object
            if (length(rownames)) 
                dimnames(ans$X)[[1]] <- rownames[ok]
            else {
                xx <- seq(1, length(na.x))
                dimnames(ans$X) <- list(NULL, NULL)
                dimnames(ans$X)[[1]] <- xx[ok]
            }
            ans$fitted.values <- rep(NA, length(na.y))
            ans$fitted.values[ok] <- fitted
            names(ans$fitted.values) <- rownames
            if (qr.out) 
                ans$qr <- z$qr
            class(ans) <- "lts"
            attr(ans, "call") <- sys.call()
            return(ans)
        }
    }
    coefs <- rep(NA, p)
    names(coefs) <- xn
    if(qr.out)
        qrx <- qr(x)
    else 
        qrx <- qr(x)[c("rank", "pivot")]

    rk <- qrx$rank
    if (rk < p) {
        stop("x is singular")
    }
    else 
        piv <- 1:p

    if (!length(alpha)) {
        alpha <- 1/2
        quan <- quan.f(alpha, n, rk)
    }else {
        if (alpha < 1/2) 
            stop("alpha is out of range!")
        quan <- quan.f(alpha, n, rk)
    }
    y <- as.matrix(y)
    x1 <- matrix(0, ncol = p + 1, nrow = n)
    x1 <- cbind(x, y)
    x1 <- as.matrix(x1)
    storage.mode(x1) <- "double"
    datt <- matrix(0, ncol = p + 1, nrow = n)
    storage.mode(datt) <- "double"
    nvad <- p + 1
    inbest <- matrix(10000, nrow = quan, ncol = 1)
    storage.mode(inbest) <- "integer"
    objfct <- 0

    interc <- ifelse(intercept, 1, 0)
    intadjust <- ifelse(adjust, 1, 0)

    storage.mode(interc) <- "integer"
    storage.mode(seed) <- "integer"
    z <- .Fortran("rfltsreg", 
                x1 = x1, 
                as.integer(n), 
                as.integer(p), 
                as.integer(quan), 
                as.integer(nsamp), 
                inbest = inbest, 
                objfct = as.double(objfct), 
                as.integer(interc), 
                as.integer(intadjust),
                as.integer(nvad), 
                datt,
                seed,
                PACKAGE="rrcov")

    # vt:: lm.fit.qr == lm.fit(...,method=qr,...)
    #  cf <- lm.fit.qr(x[z$inbest, , drop = FALSE], y[z$inbest])$coef
    cf <- lm.fit(x[z$inbest, , drop = FALSE], y[z$inbest])$coef
    
    ans$best <- sort(as.vector(z$inbest))
    fitted <- x %*% cf
    resid <- y - fitted
    coefs[piv] <- cf
    ans$raw.coefficients[2:p] <- coefs[1:(p - 1)]
    ans$raw.coefficients[1] <- coefs[p]
    names(ans$raw.coefficients)[2:p] <- names(coefs)[1:(p - 1)]
    names(ans$raw.coefficients)[1] <- names(coefs)[p]
    ans$alpha <- alpha
    ans$quan <- quan
    correct <- correctiefactor.s(p, intercept = intercept, n, 
        alpha)
    s0 <- sqrt((1/quan) * sum(sort(resid^2, quan)[1:quan]))
    sh0 <- s0
    s0 <- s0 * (1/sqrt(1 - ((2 * n)/(quan * (1/qnorm((quan + 
        n)/(2 * n))))) * dnorm(1/(1/(qnorm((quan + n)/(2 * n))))))) * 
        correct
    weights <- rep(NA, n)
    if (abs(s0) < 1e-07) {
        weights <- ifelse(abs(resid) <= 1e-07, 1, 0)
        ans$scale <- ans$raw.scale <- 0
        ans$coefficients <- ans$raw.coefficients
    }
    else {
        ans$raw.scale <- s0
        ans$raw.resid <- resid/ans$raw.scale
        quantiel <- qnorm(0.9875)
        weights <- ifelse(abs(resid/s0) <= quantiel, 1, 0)
                
        # vt:: weights has to be a vector instead of a matrix - 
        #      to avoid "Error in x * wtmult : non-conformable arrays"
        # 
        weights <- as.vector(weights)
        
        z1 <- lsfit(x, y, wt = weights, intercept = FALSE)
        ans$coefficients[2:p] <- z1$coef[1:(p - 1)]
        ans$coefficients[1] <- z1$coef[p]
        fitted <- x %*% z1$coef
        resid <- z1$residuals
        ans$scale <- sqrt(sum(weights * resid^2)/(sum(weights) - 
            1))
        if (sum(weights) == n) {
            cdelta.rew <- 1
            correct.rew <- 1
        }
        else {
            cdelta.rew <- (1/sqrt(1 - ((2 * n)/(sum(weights) * 
                (1/qnorm((sum(weights) + n)/(2 * n))))) * dnorm(1/(1/(qnorm((sum(weights) + 
                n)/(2 * n)))))))
            correct.rew <- correctiefactor.rew.s(p, intercept = intercept, 
                n, alpha)
        }
        ans$scale <- ans$scale * cdelta.rew * correct.rew
        ans$resid <- resid/ans$scale
        quantiel <- qnorm(0.9875)
        weights <- ifelse(abs(resid/ans$scale) <= quantiel, 1, 
            0)
    }
    names(ans$coefficients) <- names(ans$raw.coefficients)
    ans$lts.wt <- matrix(NA, length(na.y))
    ans$lts.wt[ok] <- weights
    ans$crit <- z$objfct
    if (intercept) {
        initcov <- matrix(0, nrow = 1, ncol = 1)
        initmean <- matrix(0, nrow = 1, ncol = 1)
        inbest <- matrix(0, nrow = quan, ncol = 1)
        weights <- matrix(0, nrow = n, ncol = 1)
        coeff <- matrix(0, nrow = 5, ncol = 1)
        adjustcov <- matrix(0, nrow = 1, ncol = 1)
        storage.mode(initcov) <- "double"
        storage.mode(initmean) <- "double"
        storage.mode(inbest) <- "integer"
        storage.mode(weights) <- "integer"
        storage.mode(coeff) <- "double"
        storage.mode(adjustcov) <- "double"
        storage.mode(y) <- "double"
        storage.mode(seed) <- "integer"
        exactfit <- 0
        k <- 1
        sh <- .Fortran("rffastmcd", 
                    as.matrix(y), 
                    as.integer(n), 
                    as.integer(k), 
                    as.integer(quan), 
                    nsamp = 0, 
                    initcovariance = initcov, 
                    initmean = initmean, 
                    inbest = inbest, 
                    mcdestimate = 0, 
                    weights = weights, 
                    as.integer(exactfit), 
                    coeff = coeff, 
                    kount = 0, 
                    adjustcov = adjustcov,
                    seed,
                    PACKAGE="rrcov")
        y <- as.vector(y)
        sh <- as.double(sh$adjustcov)
        ans$rsquared <- 1 - (sh0/sh)^2
    }
    else {
        s1 <- sum(sort(resid^2, quan)[1:quan])
        sh <- sum(sort(y^2, quan)[1:quan])
        ans$rsquared <- 1 - (s1/sh)
    }
    if (ans$rsquared > 1) {
        ans$rsquared <- 1
    }
    if (ans$rsquared < 0) {
        ans$rsquared <- 0
    }
    attributes(resid) <- attributes(fitted) <- attributes(y)
    ans$residuals <- rep(NA, length(na.y))
    ans$residuals[ok] <- resid
    ans$intercept <- intercept
    ans$method <- paste("Least Trimmed Squares Robust Regression.")
    if (abs(s0) < 1e-07) 
        ans$method <- paste(ans$method, , "\nAn exact fit was found!")
    if (mcd) {

# vt:: changed name of the function
#       mcd <- cov.mcd.default(X, alpha = alpha, print.it = FALSE)
        mcd <- covMcd(X, alpha = alpha, print.it = FALSE)
        if (-(determinant(mcd$cov, log = TRUE)$modulus[1] - 0)/p > 
            50) {
            ans$RD[1] <- "singularity"
        }
        else {
            ans$RD <- rep(NA, length(na.y))
            ans$RD[ok] <- sqrt(mahalanobis(X, mcd$center, mcd$cov))
            names(ans$RD) <- rownames
        }
    }
    names(ans$residuals) <- rownames
    names(ans$lts.wt) <- rownames
    names(ans$scale) <- names(ans$raw.scale) <- yn
    names(ans$rsquared) <- yn
    names(ans$crit) <- yn
    ans$X <- x
    ans$Y <- y          # VT:: 01.09.2004 - add y to the result object
    if (length(rownames)) 
        dimnames(ans$X)[[1]] <- rownames[ok]
    else {
        xx <- seq(1, length(na.x))
        dimnames(ans$X) <- list(NULL, NULL)
        dimnames(ans$X)[[1]] <- xx[ok]
    }
    ans$fitted.values <- rep(NA, length(na.y))
    ans$fitted.values[ok] <- fitted
    names(ans$fitted.values) <- rownames
    if (qr.out) 
        ans$qr <- qrx
    class(ans) <- "lts"
    attr(ans, "call") <- sys.call()
    return(ans)
}

#predict.lts <- function (object, newdata, na.action = na.pass, ...)
#{
#    if (missing(newdata)) return(fitted(object))
#    ## work hard to predict NA for rows with missing data
#    Terms <- delete.response(terms(object))
#    m <- model.frame(Terms, newdata, na.action = na.action,
#                     xlev = object$xlevels)
#    if(!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m)
#    X <- model.matrix(Terms, m, contrasts = object$contrasts)
#    drop(X %*% object$coefficients)
#} 
##  rrcov : Scalable Robust Estimators with High Breakdown Point
##
##  This program is free software; you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation; either version 2 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program; if not, write to the Free Software
##  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
##

print.lts <- function (x, digits = max(3, getOption("digits") - 3), ...)
{
    if(!is.null(cl <- x$call)) {
        cat("Call:\n")
        dput(cl)
        cat("\n")
    }
    cat("Coefficients:\n")
    print.default(format(coef(x), digits = digits), print.gap = 2, quote = FALSE)
    cat("\nScale estimates", format(x$scale, digits = digits) ,"\n\n")
    invisible(x)
}
##  rrcov : Scalable Robust Estimators with High Breakdown Point
##
##  This program is free software; you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation; either version 2 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program; if not, write to the Free Software
##  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
##
print.mcd <- function (x, digits = max(3, getOption("digits") - 3), ...)
{
    if(!is.null(cl <- x$call)) {
        cat("Call:\n")
        dput(cl)
        cat("\n")
    }
    cat("\nLog(det): ", format(log(x$crit), digits = digits) ,"\n\n")
    cat("Center:\n")
    print.default(format(x$center, digits = digits), print.gap = 2, quote = FALSE)
    cat("\nCovariance Matrix:\n")
    print.default(format(x$cov, digits = digits), print.gap = 2, quote = FALSE)
    invisible(x)
}
##  rrcov : Scalable Robust Estimators with High Breakdown Point
##
##  This program is free software; you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation; either version 2 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program; if not, write to the Free Software
##  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
#
#   I would like to thank Peter Filtzmoser for providing the initial code of 
#   this function.


tolellipse <- function(x, 
                        mcd, 
                        cutoff, 
                        id.n,
                        classic=FALSE,
                        tol.inv=1e-07) {

##@bdescr
## Tolerance Ellipse Plot: 
##    Plots the 97.5% tolerance ellipse of the bivariate data set (x).
##    The ellipse is defined by those data points whose distance (dist) 
##    is equal to the squareroot of the 97.5% chisquare quantile with 
##    2 degrees of freedom. 

##@edescr
##
##@in  x                 : [matrix] A data.frame or matrix, n > 2*p
##@in  mcd               : [mcd object] An object of type mcd - its attributes 
##                                      center and cov will be used
##@in  cutoff            : [number] Distance needed to flag data points outside the ellipse 
##@in  outflag           : [logical] Whether to print the labels of the outliers 
##@in  tol.inv           : [number] tolerance to be used for computing the inverse see 'solve'.
##                                  defaults to 1e-7 

    
    ellips <- function(x, y, loc, cov) {
        # calculates a 97,5% ellipsoid 
        # input: data set, location and covariance estimate, cutoff
    
        dist <- sqrt(qchisq(0.975, 2))
        A <- solve(cov)
        lambda1 <- max(eigen(A)$values)
        lambda2 <- min(eigen(A)$values)
        eigvect <- eigen(A)$vectors[, order(eigen(A)$values)[2]]
        z <- seq(0, 2 * pi, 0.01)
        z1 <- dist/sqrt(lambda1) * cos(z)
        z2 <- dist/sqrt(lambda2) * sin(z)
        alfa <- atan(eigvect[2]/eigvect[1])
        r <- matrix(c(cos(alfa),  - sin(alfa), sin(alfa), cos(alfa)), ncol = 2)
        z <- t(t(cbind(z1, z2) %*% r) + loc)    #   xmin <- min(x, z[, 1])
    
    #   xmax <- max(x, z[, 1])
    #   ymin <- min(y, z[, 2])
    #   ymax <- max(y, z[, 2])
    #   print(xmin)
    #   print(xmax)
    #   print(ymin)
    #   print(ymax)
    #   plot(x, y, xlim = c(xmin, xmax), ylim = c(ymin, ymax), type = "n")
    #   points(z[, 1], z[, 2], type = "l")
    #   points(x, y)
        z
    }

    ##  parameters and preconditions     

    if(is.vector(x) || is.matrix(x)) {
        if(!is.numeric(x))
            stop(message = "x is not a numeric dataframe or matrix.")
    }else if(is.data.frame(x)) {
        if(!all(sapply(x,data.class) == "numeric"))
            stop(message = "x is not a numeric dataframe or matrix.")
    }

    n <- dim(x)[1]
    p <- dim(x)[2]
    
    if(missing(cutoff))
        cutoff <- sqrt(qchisq(0.975, 2))

    if(p != 2)
        stop("Dimension must be 2!")
    
    if(missing(mcd))
        mcd <- covMcd(x)

    if(length(mcd$center)  == 0 ||  length(mcd$cov) == 0)
        stop(message = "Invalid mcd object: attributes center and cov missing!")

    x.loc <- mcd$center
    x.cov <- n/(n - 1) * mcd$cov
    z1 <- ellips(x[, 1], x[, 2], loc = apply(x, 2, mean), cov = n/(n - 1) * cov.wt(x)$cov)
    z2 <- ellips(x[, 1], x[, 2], loc = x.loc, cov = x.cov)
    x1 <- c(min(x[, 1], z1[, 1], z2[, 1]), max(x[,1],z1[,1], z2[,1]))
    y1 <- c(min(x[, 2], z1[, 2], z2[, 2]), max(x[,2],z1[,2], z2[,2]))
    
    md <- sqrt(mahalanobis(x,apply(x,2,mean),cov(x), tol.inv=tol.inv))
    rd <- sqrt(mahalanobis(x,mcd$center,mcd$cov, tol.inv=tol.inv))
    
    if(classic)
        par(mfrow = c(1, 2))
    else
        par(mfrow = c(1, 1))
    
    if(missing(id.n))
        id.n <- length(which(rd>cutoff))          
    ind <- sort(rd, index.return=TRUE)$ix
    ind <- ind[(n-id.n+1):n]

##  1. Robust tollerance
##  define the plot, plot a box, plot the "good" points, 
##  plot the outliers either as points or as numbers depending on outflag,
##  plot the ellipse, write a title of the plot
    plot(x[, 1], x[, 2], xlim = x1, ylim = y1, xlab = "", ylab = "", type = "p")
    box()
    xrange <- par("usr")
    xrange <- xrange[2] - xrange[1]
    text(x[ind, 1] + xrange/50, x[ind, 2], ind)

    points(z2[, 1], z2[, 2], type = "l")
    title(main = "ROBUST TOLERANCE \n    ELLIPSE (97.5%)")

##  2. Classical tollerance
    if(classic){
        plot(x[, 1], x[, 2], xlim = x1, ylim = y1, xlab = "", ylab = "", type = "p")
        box()
        
        xrange <- par("usr")
        xrange <- xrange[2] - xrange[1]
        text(x[ind, 1] + xrange/50, x[ind, 2], ind)
        
        points(z1[, 1], z1[, 2], type = "l")
        title(main = "CLASSICAL TOLERANCE \n    ELLIPSE (97.5%)")
    }
    invisible()
}
.First.lib <- function(lib, pkg) {

    where <- match(paste("package:", pkg, sep = ""), search())
    ver <- read.dcf(file.path(lib, pkg, "DESCRIPTION"), "Version")
    ver <- as.character(ver)
    title <- read.dcf(file.path(lib, pkg, "DESCRIPTION"), "Title")
    title <- as.character(title)
    cat(paste(title, " (version ", ver, ")\n", sep = ""))
    
    library.dynam("rrcov", pkg, lib)
}
