f.tau.huber <- function(dat, loc, tc=1.345){
  ## Purpose: Korrekturfaktor Tau für die Varianz des Huber-M-Schätzern
  ## -------------------------------------------------------------------------
  ## Arguments: dat =  Daten
  ##            loc =  Lokations-Punkt
  ##            tc  =  Parameter der Huber Psi-Funktion
  ## -------------------------------------------------------------------------
  ## Author: Rene Locher
  ## Update: R. Frisullo 23.4.02 (Kommentar eingeführt)
  resid <- (dat - loc)/mad(dat)
  psi <- ifelse(abs(resid) <= tc, resid, sign(resid)*tc)
  psiP <- ifelse(abs(resid) <= tc, 1, 0)
  return(length(dat) * sum(psi^2) / sum(psiP)^2)
}


p.profileTraces <- function (x, cex = 1){
  ## Purpose: Plot a profile.nls Object
  ## -------------------------------------------------------------------------
  ## Arguments: x =  profile.nls Object
  ##            
  ## Description:  Displays a series of plots of the profile t
  ##               function and the likelihood profile traces for the
  ##               parameters in a nonlinear regression model that has
  ##               been fitted with `nls' and profiled with
  ##               `profile.nls'
  ##
  ## -------------------------------------------------------------------------
  ## Author: Andreas Ruckstuhl
  ##         R port by Isabelle Flückiger and Marcel Wolbers
    nx <- names(x)
    np <- length(x)
    opar <- par(oma = c(2, 2, 1.5, 0), mfrow = c(np, np), mar = c(2, 
        4, 0, 0) + 0.2)
    on.exit(par(opar))
    for (i in 1:np) {
        for (j in 1:i) {
            if (i == j) {
                if (!is.null(this.comp <- x[[i]])) {
                  xx <- this.comp$par[, nx[i]]
                  tau <- this.comp[[1]]
                  plot(spline(xx, tau), xlab = "", ylab = "", 
                    type = "l", las = 1, mgp = c(3, 0.8, 0), 
                    cex = 0.5 * cex)
                  points(xx[tau == 0], 0, pch = 3)
                  pusr <- par("usr")
                  if (is.R()) {
                    mtext(side = 1, line = 0.8, at = -1/(2 * 
                      np) + i/np, text = nx[j], outer = TRUE, 
                      cex = cex)
                    mtext(side = 2, line = 0.8, at = 1 + 1/(2 * 
                      np) - i/np, text = nx[i], outer = TRUE, 
                      cex = cex)
                  }
                  else {
                    mtext(side = 1, line = 0.8, at = mean(pusr[1:2]), 
                      text = nx[j], outer = TRUE, cex = cex)
                    mtext(side = 2, line = 0.8, at = mean(pusr[3:4]), 
                      text = nx[i], outer = TRUE, cex = cex)
                  }
                }
            }
            else {
                if ((!is.null(x.comp <- x[[j]])) & (!is.null(y.comp <- x[[i]]))) {
                  xx <- x.comp$par[, nx[j]]
                  xy <- x.comp$par[, nx[i]]
                  yx <- y.comp$par[, nx[j]]
                  yy <- y.comp$par[, nx[i]]
                  plot(xx, xy, xlab = "", ylab = "", las = 1, 
                    mgp = c(3, 0.8, 0), type = "n", xlim = range(c(xx, 
                      yx)), ylim = range(c(xy, yy)), cex = 0.5 * 
                      cex)
                  lines(xx, xy, col = 2)
                  lines(yx, yy, col = 3)
                }
            }
        }
        if (i < np) 
            for (k in 1:(np - i + if (is.R()) 
                0
            else 1)) frame()
    }
    mtext(side = 3, line = 0.2, text = "t-Profil-Plot und Profilspuren", 
        outer = TRUE, cex = 1.2 * cex)
}


p.Serie2Aufg2c <- function (x, cex = 1){
  ## Purpose: Plot a profile.nls Object
  ## -------------------------------------------------------------------------
  ## Arguments: x =  profile.nls Object
  ##
  ## Description:  Modification of the R-function p.profileTraces:
  ##               Displays a series of plots of the profile t
  ##               function and the likelihood profile traces for the 
  ##               parameters - without the first -
  ##               in a nonlinear regression model that has
  ##               been fitted with `nls' and profiled with
  ##               `profile.nls'
  ##
  ## -------------------------------------------------------------------------
  ## Author of p.profileTraces: Andreas Ruckstuhl
  ##                            modified by Roberto Frisullo
  ##         
    nx <- names(x)
    np <- length(x) 
    opar <- par(oma = c(2, 2, 1.5, 0), mfrow = c(np-1, np-1), mar = c(2, 
        4, 0, 0) + 0.2)
    on.exit(par(opar))
    for (i in 2:np) {
        for (j in 2:i) {
            if (i == j) {
                if (!is.null(this.comp <- x[[i]])) {
                  xx <- this.comp$par[, nx[i]]
                  tau <- this.comp[[1]]
                  plot(spline(xx, tau), xlab = "", ylab = "", 
                    type = "l", las = 1, mgp = c(3, 0.8, 0), 
                    cex = 0.5 * cex)
                  points(xx[tau == 0], 0, pch = 3)
                  pusr <- par("usr")
                  if (is.R()) {
                    mtext(side = 1, line = 0.8, at = -1/(2 * 
                      (np-1)) + (i-1)/(np-1), text = nx[j], outer = TRUE, 
                      cex = cex)
                    mtext(side = 2, line = 0.8, at = 1 + 1/(2 * 
                      (np-1)) - (i-1)/(np-1), text = nx[i], outer = TRUE, 
                      cex = cex)
                  }
                  else {
                    mtext(side = 1, line = 0.8, at = mean(pusr[1:2]), 
                      text = nx[j], outer = TRUE, cex = cex)
                    mtext(side = 2, line = 0.8, at = mean(pusr[3:4]), 
                      text = nx[i], outer = TRUE, cex = cex)
                  }
                }
            }
            else {
                if ((!is.null(x.comp <- x[[j]])) & (!is.null(y.comp <- x[[i]]))) {
                  xx <- x.comp$par[, nx[j]]
                  xy <- x.comp$par[, nx[i]]
                  yx <- y.comp$par[, nx[j]]
                  yy <- y.comp$par[, nx[i]]
                  plot(xx, xy, xlab = "", ylab = "", las = 1, 
                    mgp = c(3, 0.8, 0), type = "n", xlim = range(c(xx, 
                      yx)), ylim = range(c(xy, yy)), cex = 0.5 * 
                      cex)
                  lines(xx, xy, col = 2)
                  lines(yx, yy, col = 3)
                }
            }
        }
        if (i < np) 
            for (k in 1:(np - i + if (is.R()) 
                0
            else 1)) frame()
    }
    mtext(side = 3, line = 0.2, text = "t-Profil-Plot und Profilspuren", 
        outer = TRUE, cex = 1.2 * cex)
}


f.calib <- function(cal.x, cal.pred, cal.limit, cal.y){
  ## Purpose: Bestimmung des Kalibrationsintervalls
  ## -------------------------------------------------------------------------
  ## Arguments: cal.x:      x-Werte der Kalibrationskurve                      
  ##            cal.pred    predict-Objekt fuer Kalibrationskurve              
  ##            cal.limit   Halbe Breite des Vorhersageintervalls              
  ##            cal.y:      beobachteter y-Wert, zu dem der x-Wert gesucht wird
  ##
  ## -------------------------------------------------------------------------
  ## Author: Andreas Ruckstuhl, Roberto Frisullo

  I1 <- uniroot(function(x) cal.y - 
                approx(cal.x,cal.pred-cal.limit,xout=x)$y,
  		       range(cal.x))$root
  Ic <- uniroot(function(x) cal.y - 
                approx(cal.x,cal.pred,xout=x)$y,
  		       range(cal.x))$root
  I2 <- uniroot(function(x) cal.y -
           approx(cal.x,cal.pred+cal.limit,xout=x)$y,
  		       range(cal.x))$root

   data.frame("lower"=min(I1,I2),"center"=Ic,"upper"=max(I1,I2))
}


p.ldv <- function(ldafit, data=D[,2:3], group=D$Klasse){
  ## Rks, Juni 02
  ldaDV <- predict(ldafit, data, dimen=2)$x
  eqscplot(ldaDV[,1], ldaDV[,2], type="n",
	 xlab="first discriminant variable",
	 ylab="second discriminant variable", main="", las=1)
  for(i in 1:4){
    x.c <- c("A", "B", "C", "D")[i]; ok <- group==x.c
    text(ldaDV[ok,1], ldaDV[ok,2], labels=x.c,
         cex=0.6, col=i+1)
  }
  ## Add centre of groups
  x <- apply(diag(ldafit$prior) %*% ldafit$means, 2, sum) 
  z <- scale(ldafit$means, x, F)      ## remove overall means
  points(z %*% ldafit$scaling , col=2:5, pch=20)
  points(z %*% ldafit$scaling, pch=1)
}


p.predplot <- function(ldafit, data=D[,2:4], group=D$Klasse, len=100){
  ## Rks, Juni 02
  ## From V&R, 2nd ed
  g <- as.factor(group)
  lev <- levels(g)
  plot(data[,1], data[,2], type="n", las=1)
  text(data[,1], data[,2], labels=as.character(g), col=codes(g)+1, cex=0.6)
  h <- range(data[,1]); xp <- seq(h[1], h[2], length=len)
  h <- range(data[,2]); yp <- seq(h[1], h[2], length=len)
  xyp <- expand.grid(x1=xp, x2=yp)
  Z <- predict(ldafit, newdata=xyp)
  zp <- Z$post[,2] - pmax(Z$post[,1],Z$post[,3],Z$post[,4])
  contour(xp, yp, matrix(zp, len), add=T, levels=0, drawlabels=F)
  zp <- Z$post[,3] - pmax(Z$post[,1],Z$post[,2],Z$post[,4])
  contour(xp, yp, matrix(zp, len), add=T, levels=0, drawlabels=F)
  zp <- Z$post[,4] - pmax(Z$post[,1],Z$post[,2],Z$post[,3])
  contour(xp, yp, matrix(zp, len), add=T, levels=0, drawlabels=F)
  invisible()
}

rlda <- function(x, grouping, prior=proportions,
                         tol=1e-04, nu=5, ...){
  ## Rks, Juni 02
  ## modified version of lda.default: Estimation with cov.mcd and
  ## robust estimation of location
    if (is.null(dim(x))) 
        stop("x is not a matrix")
    x <- as.matrix(x)
    n <- nrow(x)
    p <- ncol(x)
    if (n != length(grouping)) 
        stop("nrow(x) and length(grouping) are different")
    g <- as.factor(grouping)
    lev <- lev1 <- levels(g)
    counts <- as.vector(table(g))
    if (any(counts == 0)) {
        warning(paste("group(s)", paste(lev[counts == 0], collapse = " "), 
            "are empty"))
        lev1 <- lev[counts > 0]
        g <- factor(g, levels = lev1)
        counts <- as.vector(table(g))
    }
    proportions <- counts/n
    ng <- length(proportions)
    if (any(prior < 0) || round(sum(prior), 5) != 1) 
        stop("invalid prior")
    if (length(prior) != ng) 
        stop("prior is of incorrect length")
    names(prior) <- names(counts) <- lev1
     ## robust estiamtion of group location
    group.means <- matrix(NA, ncol=ncol(x), nrow=length(lev1),
                          dimnames=list(lev1,names(x))) 
    for(i in lev1)
      group.means[i,] <- cov.rob(x[g==i,], method="mcd")$center
    f1 <- sqrt(diag(var(x - group.means[g, ])))
    if (any(f1 < tol)) 
        stop(paste("variable(s)", paste(format((1:p)[f1 < tol]), 
            collapse = " "), "appear to be constant within groups"))
    scaling <- diag(1/f1, , p)
    cov <- n/(n - ng) * cov.rob((x - group.means[g, ]) %*% scaling,
                                method="mcd")$cov
    sX <- svd(cov, nu = 0)
    rank <- sum(sX$d > tol^2)
    if (rank < p) 
      warning("variables are collinear")
    scaling <- scaling %*% sX$v[, 1:rank] %*% diag(sqrt(1/sX$d[1:rank]), 
                                                   , rank)
    ##
    xbar <- apply(prior %*% group.means, 2, sum)
    fac <- 1/(ng - 1)
    X <- sqrt((n * prior) * fac) * scale(group.means, center = xbar, 
        scale = FALSE) %*% scaling
    X.s <- svd(X, nu = 0)
    rank <- sum(X.s$d > tol * X.s$d[1])
    scaling <- scaling %*% X.s$v[, 1:rank]
    if (is.null(dimnames(x))) 
        dimnames(scaling) <- list(NULL, paste("LD", 1:rank, sep = ""))
    else {
        dimnames(scaling) <- list(colnames(x), paste("LD", 1:rank, 
            sep = ""))
        dimnames(group.means)[[2]] <- colnames(x)
    }
    structure(list(prior = prior, counts = counts, means = group.means, 
                   scaling = scaling, lev = lev, svd = X.s$d[1:rank], N = n, 
                   call = match.call()), class = "lda")
}
