#   translated biweight s-estimation 
#   corrected (?) version
#   corrections made in functions rho.bt, psi.bt, wt.bt, erho.bt
#   functions added are erho.bt.lim, erho.bt.lim.p
#
rho.bt <- function(x,c1,M)
{
    if (c1 == 0)
    {
    x1 <- x - M
    ivec1 <- (x1 <= 0)
    ivec2 <- (x1 > 0)
    }
    else
    {
    x1 <- (x-M)/c1
    ivec1 <- (x1 < 0)
    ivec2 <- (x1 >  1)
    }
    return(ivec1*(x^2/2)
        +ivec2*(M^2/2+c1*(5*c1+16*M)/30)
        +(1-ivec1-ivec2)*(M^2/2-M^2*(M^4-5*M^2*c1^2+15*c1^4)/(30*c1^4)
            +(1/2+M^4/(2*c1^4)-M^2/c1^2)*x^2
            +(4*M/(3*c1^2)-4*M^3/(3*c1^4))*x^3
            +(3*M^2/(2*c1^4)-1/(2*c1^2))*x^4
            -4*M*x^5/(5*c1^4)+x^6/(6*c1^4)))
}
psi.bt <- function(x,c1,M)
{
    if (c1 == 0)
    {
    x1 <- x - M
    ivec1 <- (x1 <= 0)
    ivec2 <- (x1 > 0)
    }
    else
    {
    x1 <- (x-M)/c1
    ivec1 <- (x1 < 0)
    ivec2 <- (x1 >  1)
    }
    return(ivec1*x+(1-ivec1-ivec2)*x*(1-x1^2)^2)
}
psip.bt <- function(x,c1,M)
{
    x1 <- (x-M)/c1
    ivec1 <- (x1 < 0)
    ivec2 <- (x1 >  1)
    return(ivec1+(1-ivec1-ivec2)*((1-x1^2)^2-4*x*x1*(1-x1^2)/c1))
}
wt.bt <- function(x,c1,M)
{
    if (c1 == 0)
    {
    x1 <- x - M
    ivec1 <- (x1 <= 0)
    ivec2 <- (x1 > 0)
    }
    else
    {
    x1 <- (x-M)/c1
    ivec1 <- (x1 < 0)
    ivec2 <- (x1 >  1)
    }
    return(ivec1+(1-ivec1-ivec2)*(1-x1^2)^2)
}
v.bt <- function(x,c1,M) return(x*psi.bt(x,c1,M)) 
#
#       functions for computing expectations
#
chi.int <- function(p,a,c1)	#
#
#   partial expectation d in (0,c1) of d^a under chi-squared p	#
#
  return( exp(lgamma((p+a)/2)-lgamma(p/2))*2^{a/2}*pchisq(c1^2,p+a) )
#
#
chi.int2 <- function(p,a,c1)	#
#
#   partial expectation d in (c1,\infty) of d^a under chi-squared p	#
#
  return( exp(lgamma((p+a)/2)-lgamma(p/2))*2^{a/2}*(1-pchisq(c1^2,p+a)) )
#
#
#       derivatives of above wrt c1
#
chi.int.p <- function(p,a,c1)
  return( exp(lgamma((p+a)/2)-lgamma(p/2))*2^{a/2}*dchisq(c1^2,p+a)*2*c1 )
#
#
chi.int2.p <- function(p,a,c1)
  return( -exp(lgamma((p+a)/2)-lgamma(p/2))*2^{a/2}*dchisq(c1^2,p+a)*2*c1 )
#
#       functions for expectations with translated biweight
#
epsi.bt <- function(p,c1,M)	#
#
#   expectation of psi(d) under chi-squared p	#
#
  return(chi.int(p,1,M)
        +(1-2*M^2/c1^2+M^4/c1^4)*(chi.int(p,1,M+c1)-chi.int(p,1,M))
        +(4*M/c1^2-4*M^3/c1^4)*(chi.int(p,2,M+c1)-chi.int(p,2,M))
        +(6*M^2/c1^4-2/c1^2)*(chi.int(p,3,M+c1)-chi.int(p,3,M))
        -(4*M/c1^4)*(chi.int(p,4,M+c1)-chi.int(p,4,M))
        +(chi.int(p,5,M+c1)-chi.int(p,5,M))/c1^4)	
#
#
epsip.bt <- function(p,c1,M)	#
#
#   expectation of derivative of psi(d) under chi-squared p	#
#
  return(chi.int(p,0,M)
        +(1-2*M^2/c1^2+M^4/c1^4)*(chi.int(p,0,M+c1)-chi.int(p,0,M))
        +(4*M/c1^2-4*M^3/c1^4)*(chi.int(p,1,M+c1)-chi.int(p,1,M))*2
        +(6*M^2/c1^4-2/c1^2)*(chi.int(p,2,M+c1)-chi.int(p,2,M))*3
        -(4*M/c1^4)*(chi.int(p,3,M+c1)-chi.int(p,3,M))*4
        +(chi.int(p,4,M+c1)-chi.int(p,4,M))*5/c1^4)	#
#
#
epsi2.bt <- function(p,c1,M)	#
#
#   expectation of psi^2(d) under chi-squared p	#
#
  return(chi.int(p,2,M)
        +(1-4*M^2/c1^2+6*M^4/c1^4-4*M^6/c1^6+M^8/c1^8)*(chi.int(p,2,M+c1)-chi.int(p,2,M))
        +(8*M/c1^2-24*M^3/c1^4+24*M^5/c1^6-8*M^7/c1^8)*(chi.int(p,3,M+c1)-chi.int(p,3,M))
        +(-4/c1^2+36*M^2/c1^4-60*M^4/c1^6+28*M^6/c1^8)*(chi.int(p,4,M+c1)-chi.int(p,4,M))
        +(-24*M/c1^4+80*M^3/c1^6-56*M^5/c1^8)*(chi.int(p,5,M+c1)-chi.int(p,5,M))
        +(6/c1^4-60*M^2/c1^6+70*M^4/c1^8)*(chi.int(p,6,M+c1)-chi.int(p,6,M))
        +(24*M/c1^6-56*M^3/c1^8)*(chi.int(p,7,M+c1)-chi.int(p,7,M))
        +(-4/c1^6+28*M^2/c1^8)*(chi.int(p,8,M+c1)-chi.int(p,8,M))
        -(8*M/c1^8)*(chi.int(p,9,M+c1)-chi.int(p,9,M))
        +(1/c1^8)*(chi.int(p,10,M+c1)-chi.int(p,10,M))
        )
#
#
ew.bt <- function(p,c1,M)	#
#
#   expectation of w(d)=psi(d)/d under chi-squared p
#
  return(chi.int(p,0,M)
        +(1-2*M^2/c1^2+M^4/c1^4)*(chi.int(p,0,M+c1)-chi.int(p,0,M))
        +(4*M/c1^2-4*M^3/c1^4)*(chi.int(p,1,M+c1)-chi.int(p,1,M))
        +(6*M^2/c1^4-2/c1^2)*(chi.int(p,2,M+c1)-chi.int(p,2,M))
        -(4*M/c1^4)*(chi.int(p,3,M+c1)-chi.int(p,3,M))
        +(chi.int(p,4,M+c1)-chi.int(p,4,M))/c1^4)	#
#
#
erho.bt <- function(p,c1,M)	#
#
#   expectation of rho(d) under chi-squared p
#
{
    if (c1 == 0)
    {
    return(chi.int(p,2,M)/2 + M^2/2*chi.int2(p,0,M))
    }
    else
    {
    return(chi.int(p,2,M)/2
        +(M^2/2+c1*(5*c1+16*M)/30)*chi.int2(p,0,M+c1)
        +(M^2/2-M^2*(M^4-5*M^2*c1^2+15*c1^4)/(30*c1^4))*(chi.int(p,0,M+c1)-chi.int(p,0,M))
        +(1/2+M^4/(2*c1^4)-M^2/c1^2)*(chi.int(p,2,M+c1)-chi.int(p,2,M))
        +(4*M/(3*c1^2)-4*M^3/(3*c1^4))*(chi.int(p,3,M+c1)-chi.int(p,3,M))
        +(3*M^2/(2*c1^4)-1/(2*c1^2))*(chi.int(p,4,M+c1)-chi.int(p,4,M))
        -(4*M/(5*c1^4))*(chi.int(p,5,M+c1)-chi.int(p,5,M))
        +(1/(6*c1^4))*(chi.int(p,6,M+c1)-chi.int(p,6,M)))
    }
}
#
#
#
erho.bt.lim <- function(p,c1)
    return(chi.int(p,2,c1) + c1^2*chi.int2(p,0,c1))
#
#
erho.bt.lim.p <- function(p,c1)
    return(chi.int.p(p,2,c1) + 2*c1*chi.int2(p,0,c1)
           + c1^2*chi.int2.p(p,0,c1))
#
#
#       find constants needed for s-estimation
#
rejpt.bt.lim <- function(p,r)	#
#
#   find p-value of translated biweight limit c 
#   that gives a specified breakdown
#
{
    c1 <- 2*p
    iter <- 1
    crit <- 100
    eps <- 1e-5
    while ((crit > eps)&(iter<100))
    {
        c1.old <- c1
        fc <- erho.bt.lim(p,c1) - c1^2*r
        fcp <- erho.bt.lim.p(p,c1) - 2*c1*r
        c1 <- c1 - fc/fcp
        if (c1 < 0)  c1 <- c1.old/2
        crit <- abs(fc)
        iter <- iter+1
    }
    return(c(c1,pchisq(c1^2,p),log10(1-pchisq(c1^2,p))))
}	#
#
#
csolve.bt <- function(n,p,r,alpha,asymp=FALSE)	#
#
#   find constants c1 and M that gives a specified breakdown r
#   and rejection point alpha
#
{
    if (asymp == FALSE){if (r > (n-p)/(2*n) ) r <- (n-p)/(2*n)} # 
# maximum achievable breakdown
#
#   if rejection is not achievable, use c1=0 and best rejection
#
    limvec <- rejpt.bt.lim(p,r)
    if (1-limvec[2] <= alpha)
    {
        c1 <- 0
        M <- sqrt(qchisq(1-alpha,p))
        print("adjusting alpha")
        print(c(alpha,M,c1))
    }
    else 
    {
    c1.plus.M <- sqrt(qchisq(1-alpha,p))
    M <- sqrt(p)
    c1 <- c1.plus.M - M
    iter <- 1
    crit <- 100
    eps <- 1e-5
    while ((crit > eps)&(iter<100))
    {
        deps <- 1e-4
        M.old <- M
        c1.old <- c1
        er <- erho.bt(p,c1,M)
        fc <- er - r*(M^2/2+c1*(5*c1+16*M)/30)
        fcc1 <- (erho.bt(p,c1+deps,M)-er)/deps
        fcM  <- (erho.bt(p,c1,M+deps)-er)/deps
        fcp <- fcM - fcc1 - r*(M-(5*c1+16*M)/30+c1*11/30)
        M <- M - fc/fcp
        if (M >= c1.plus.M ){M <- (M.old + c1.plus.M)/2}
        c1 <- c1.plus.M - M
        crit <- abs(fc)
        iter <- iter+1
    }
    }
    return(list(c1=c1,M=M))
}
#
#
ksolve.bt <- function(d,p,c1,M,b0)	#
#
#       find a constant k which satisfies the s-estimation constraint
#       for modified biweight	#
#
{
    k <- 1
    iter <- 1
    crit <- 100
    eps <- 1e-5
    while ((crit > eps)&(iter<100))
    {
        k.old <- k
        fk <- mean(rho.bt(d/k,c1,M))-b0
        fkp <- -mean(psi.bt(d/k,c1,M)*d/k^2)
        k <- k - fk/fkp
        if (k < k.old/2)  k <- k.old/2
        if (k > k.old*1.5) k <- k.old*1.5
        crit <- abs(fk)
        iter <- iter+1
    }
    return(k)
}
#
#
multmest2.bt <- function(x,t1,s,c1,M,eps=1e-3,maxiter=20)	#
#
#       s-estimation from starting point (t1,s)
#       version for new s, which has no mahalanobis function
#       for modified biweight
#       ohne nicht benoetigte Berechnungen
#
{
    n <- nrow(x)
    p <- ncol(x)
    b0 <- erho.bt(p,c1,M)
    crit <- 100
    iter <- 1
    w1d <- rep(1,n)
    w2d <- w1d
    while ((crit > eps)&(iter <= maxiter))
    {
        t.old <- t1
        s.old <- s
        wt.old <- w1d
        v.old <- w2d
        d2 <- mahalanobis(x,t1,s)
        d <- sqrt(d2)
        k <- ksolve.bt(d,p,c1,M,b0)
        d <- d/k
        w1d <- wt.bt(d,c1,M)
        w2d <- v.bt(d,c1,M)
        t1 <- (w1d %*% x)/sum(w1d)
        s <- s*0
        for (i in 1:n)
        {
            xc <- as.vector(x[i,]-t1)
            s <- s + as.numeric(w1d[i])*(xc %o% xc)
        }
        s <- p*s/sum(w2d)
        crit <- max(abs(w1d-wt.old))/max(w1d)
        iter <- iter+1
    }
    return(list(t1=t1,s=s))
}
#
#
"mod.s.schaetzer"<-
function(data)
{
	N <- nrow(data)	# Anzahl der Beobachtungen
	p <- ncol(data)	# Anzahl der Variablen
	p1 <- p + 1		#
############################################################################
# MVE-Schaetzer als Startschaetzer; standard random sampling algorithm unter
# Benutzung der in S-Plus implementierten Routine fuer den MVE
###########################################################################
#
	startschaetzer <- cov.mve(data, print = F, popsize = 1, maxslen = 
		p1, births.n = 0)
	startlage <- startschaetzer[[2]]	# MVE-Lageschaetzer
	startkov <- startschaetzer[[1]]	# MVE-Kovarianzschaetzer
############################################################################
# small sample correction
############################################################################
#
	korrektur <- ((1 + 15/(N - p))^2)
	startkov <- startkov * korrektur	# korrigierter Kovarianzschaetzer
###################################################################################
# Bestimmung der Konstanten c und M fuer die Biweight-Funktion mit maximalem
# Bruchpunkt und ARP von 0.01 (nach Vorschlag aus Rocke-Artikel)
###################################################################################
#
	cons <- csolve.bt(N, p, 0.5, 0.01)
	ccons <- cons[[1]]
	Mcons <- cons[[2]]	#
###############################################################################
# Bestimmung der S-Schaetzer aus den Daten unter Verwendung der Startschaetzer
###############################################################################
#
	schaetzer <- multmest2.bt(data, startlage, startkov, ccons, Mcons)
	return(schaetzer)
}