einfsMM <- function(Y, MMloc, MMcov) {

# empirical influences for MM-estimates + PCA estimates

# --------------------------------------------------------------------

rhobiweight <- function(x,c)
{
# Computes Tukey's biweight rho function with constant c for all values in x

hulp <- x^2/2 - x^4/(2*c^2) + x^6/(6*c^4)
rho <- hulp*(abs(x)<c) + c^2/6*(abs(x)>=c)

return(rho)
}

# --------------------------------------------------------------------

rhobiweightder1 <- function(x,c)
{
# Computes Tukey's biweight psi function with constant c for all values in x

hulp <- x - 2*x^3/(c^2) + x^5/(c^4)
rho <- hulp*(abs(x)<c)

return(rho)
}

# --------------------------------------------------------------------

rhobiweightder2 <- function(x,c)
{
# Computes derivative of Tukey's biweight psi function with constant c for all values in x

hulp <- 1 - 6*x^2/(c^2) + 5*x^4/(c^4)
rho <- hulp*(abs(x)<c)

return(rho)
}

# --------------------------------------------------------------------

vecop <- function(mat) {
# performs vec-operation (stacks colums of a matrix into column-vector)

nr <- nrow(mat)
nc <- ncol(mat)

vecmat <- rep(0,nr*nc)
for (col in 1:nc) {
    startindex <- (col-1)*nr+1
    vecmat[startindex:(startindex+nr-1)] <- mat[,col]
}
return(vecmat)
}

#------------------------------------------------------------------------------#
#                function needed to determine constant in biweight             #
#------------------------------------------------------------------------------#
# (taken from Claudia Becker's Sstart0 program)

"chi.int" <- function(p, a, c1)
return(exp(lgamma((p + a)       
        #   partial expectation d in (0,c1) of d^a under chi-squared p
  /2) - lgamma(p/2)) * 2^{a/2} * pchisq(c1^2, p + a))

"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" <- function(p, a, c1)
return(exp(lgamma((p + a)       
        #   partial expectation d in (c1,\infty) of d^a under chi-squared p
  /2) - lgamma(p/2)) * 2^{a/2} * (1 - pchisq(c1^2, p + a)))

"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)

# -------------------------------------------------------------------

"csolve.bw.asymp" <- function(p, r)
{
#   find biweight c that gives a specified breakdown
        c1 <- 9
        iter <- 1
        crit <- 100
        eps <- 0.00001
        while((crit > eps) & (iter < 100)) {
                c1.old <- c1
                fc <- erho.bw(p, c1) - (c1^2 * r)/6
                fcp <- erho.bw.p(p, c1) - (c1 * r)/3
                c1 <- c1 - fc/fcp
                if(c1 < 0)
                        c1 <- c1.old/2
                crit <- abs(fc) 
                iter <- iter + 1
        }
        return(c1)
}

"erho.bw" <- function(p, c1)
  return(chi.int(p,       #   expectation of rho(d) under chi-squared p
         2, c1)/2 - chi.int(p, 4, c1)/(2 * c1^2) + chi.int(p, 6, c1)/(6 * c1^4) + (
        c1^2 * chi.int2(p, 0, c1))/6)

"erho.bw.p" <- function(p, c1)
  return(chi.int.p(p,     #   derivative of erho.bw wrt c1
    2, c1)/2 - chi.int.p(p, 4, c1)/(2 * c1^2) + (2 * chi.int(p, 4, c1))/(2 * 
        c1^3) + chi.int.p(p, 6, c1)/(6 * c1^4) - (4 * chi.int(p, 6, c1))/(
        6 * c1^5) + (c1^2 * chi.int2.p(p, 0, c1))/6 + (2 * c1 * chi.int2(p,
        0, c1))/6)

# --------------------------------------------------------------------

"sigma1.bw" <- function(p, c1)
{  
# called by csolve.bw.MM()

gamma1.1 <- chi.int(p,2,c1) - 6*chi.int(p,4,c1)/(c1^2) + 5*chi.int(p,6,c1)/(c1^4)   
gamma1.2 <- chi.int(p,2,c1) - 2*chi.int(p,4,c1)/(c1^2) + chi.int(p,6,c1)/(c1^4)
gamma1 <- ( gamma1.1 + (p+1)*gamma1.2 ) / (p+2)

sigma1.0 <- chi.int(p,4,c1) - 4*chi.int(p,6,c1)/(c1^2) + 6*chi.int(p,8,c1)/(c1^4) - 4*chi.int(p,10,c1)/(c1^6) + chi.int(p,12,c1)/(c1^8)
return( sigma1.0 / (gamma1^2) * p/(p+2) )

}

# --------------------------------------------------------------------

csolve.bw.MM <- function(p, eff)
{
# constant for second Tukey Biweight rho-function for MM, for fixed shape-efficiency 

maxit <- 1000
eps <- 10^(-8)
diff <- 10^6
#ctest <- csolve.bw.asymp(p,.5)
ctest <- -.4024 + 2.2539 * sqrt(p) # very precise approximation for c corresponding to 50% bdp
iter <- 1
while ((diff>eps) & (iter<maxit)) {
    cold <- ctest
    ctest <- cold * eff * sigma1.bw(p,cold)
    diff <- abs(cold-ctest)
    iter <- iter+1
}
return(ctest)

}

#--------------------------------------------------------------------------#
#                      end 'constant determining'                          #
#--------------------------------------------------------------------------#

# --------------------------------------------------------------------
# -                         main function                            -
# --------------------------------------------------------------------

n <- nrow(Y)
p <- ncol(Y)

bdp <- .5
c0 <- csolve.bw.asymp(p,bdp)
b <- erho.bw(p, c0)
c1 <- csolve.bw.MM(p,.95)

MMGamma <- det(MMcov)^(-1/p) * MMcov

# compute eigenvalues/vectors
eigresGamma <- eigen(MMGamma)
IXGamma <- order(eigresGamma$values)
eigs <- eigresGamma$values[IXGamma]
eigvecs <- eigresGamma$vectors[,IXGamma]

for (k in 1:p) {
    IX  <- order(abs(eigvecs[,k]))
    eigvecs[,k] <- sign(eigvecs[IX[p],k]) * eigvecs[,k]
}    

varperc <- rep(0,p-1)
for (k in 1:(p-1)) {
    varperc[k] <- sum(eigs[(p-k+1):p]) / sum(eigs)
}

divec <- sqrt(mahalanobis(Y, MMloc, MMcov))

psidervec <- rhobiweightder2(divec, c1)
psidervecS <- rhobiweightder2(divec, c0)
uvec <- rhobiweightder1(divec, c1) / divec
uvecS <- rhobiweightder1(divec, c0) / divec
vvec <- rhobiweightder1(divec,c1) * divec
vvecS <- rhobiweightder1(divec,c0) * divec
rhovecS <- rhobiweight(divec,c0)

betaMM <- (1-1/p) * mean(uvec) + 1/p * mean(psidervec)
einfs.loc <- 1 / betaMM * matrix(rep(uvec,p), ncol=p) * (Y - matrix(rep(MMloc,n), n, byrow=T))

betaS <- (1-1/p) * mean(uvecS) + 1/p * mean(psidervecS)
einfs.locS <- 1 / betaS * matrix(rep(uvecS,p), ncol=p) * (Y - matrix(rep(MMloc,n), n, byrow=T))

gamma3 <- mean(vvecS)
gamma1MM <- mean(psidervec * (divec^2) + (p+1) * vvec) / (p+2)
gamma1S <- mean(psidervecS * (divec^2) + (p+1) * vvecS) / (p+2)

einfscov <- matrix(0,n,p*p)
einfscovS <- matrix(0,n,p*p)
einfsshape <- matrix(0,n,p*p)
einfseigvec <- matrix(0,n,p*p)
einfseigscov <- matrix(0,n,p)
einfseigs <- matrix(0,n,p)
einfsvarperc <- matrix(0,n,p-1)
for (i in 1:n) {
    IFcov <- 2/gamma3 * (rhovecS[i] - b) * MMcov + 1/gamma1MM * p * vvec[i] * (t(Y[i,]-MMloc) %*% (Y[i,]-MMloc) / divec[i]^2 - 1/p * MMcov)    
    einfscov[i,] <- vecop( IFcov )
    
    IFcovS <- 2/gamma3 * (rhovecS[i] - b) * MMcov + 1/gamma1S * p * vvecS[i] * (t(Y[i,]-MMloc) %*% (Y[i,]-MMloc) / divec[i]^2 - 1/p * MMcov)
    einfscovS[i,] <- vecop( IFcovS )
    
    IFshape <- 1/gamma1MM * p * vvec[i] * det(MMcov)^(-1/p) * (t(Y[i,]-MMloc) %*% (Y[i,]-MMloc) / divec[i]^2 - 1/p * MMcov)
    einfsshape[i,] <- vecop( IFshape )
    
    IFeigvecs <- matrix(0,p,p)
    for (k in 1:p) { 
        IFeigveck <- rep(0,p)
        for (j in 1:p) {
            if (j != k) {
                IFeigveck <- IFeigveck + 1 / (eigs[k]-eigs[j]) * (t(eigvecs[,j]) %*% IFshape %*% eigvecs[,k]) %*% eigvecs[,j]
            }
        }
        IFeigvecs[,k] <- IFeigveck
        einfseigs[i,k] <- t(eigvecs[,k]) %*% IFshape %*% eigvecs[,k]
        einfseigscov[i,k] <- t(eigvecs[,k]) %*% IFcov %*% eigvecs[,k]
    }
    einfseigvec[i,] <- vecop(IFeigvecs)
    for (k in 1:(p-1)) {
        einfsvarperc[i,k] <- 1 / sum(eigs) * ((1-varperc[k]) * sum(einfseigs[i,(p-k+1):p]) - varperc[k] * sum(einfseigs[i,1:(p-k)]))
    }
}
return(list(loc=einfs.loc, locS=einfs.locS, shape=einfsshape, cov=einfscov, covS=einfscovS, eigvec=einfseigvec, eigs=einfseigs, eigscov=einfseigscov, varperc=einfsvarperc))
 
}