# R/MVTMLEr2.R In fastM: Fast Computation of Multivariate M-Estimators

```MVTMLE_symm <- function(X,nu=0,delta=10^(-7),
prewhitened=FALSE,
steps=TRUE,nmax=500)
# Computes symmetrized M-estimator of scatter.
{
Xnames <- colnames(X)
Xs <- as.matrix(X)
n <- dim(Xs)[1]
p <- dim(Xs)[2]
if (n <= p)
{
return("Too few observations!")
}
if (!prewhitened)
{
# Prewhiten the data:
S0 <- 2*cov(Xs)
B <- SqrtS(S0)
Xs <- t(qr.solve(B,t(Xs)))
}
else
{
B <- diag(rep(1,p))
Xs <- X
}

# Phase 1: Use a random subset of n
# observation pairs to get a good
# starting value:
if (steps)
{
print("Phase 1 (using a random subset of n pairs):",quote=FALSE)
}
#Pi <- sample(n)
X0 <- Xs[,] - Xs[c(2:n,1),]
C <- MVTMLE0r(X0,nu,delta,
steps=steps,prewhitened=TRUE)\$B
B <- B %*% C
Xs <- t(qr.solve(C,t(Xs)))

# Phase 2:
if (steps)
{
print("Phase 2 (using all N pairs):",quote=FALSE)
}
if (n <= nmax)
# Use a big data matrix with
# N observation pairs:
{
IJ <- combn(n,2)
XXs <- Xs[IJ[1,],] - Xs[IJ[2,],]
tmp <- MVTMLE0r(XXs,nu,delta=delta,
prewhitened=TRUE,steps=steps)
B <- B %*% tmp\$B
iter <- tmp\$iter
}
else
# Avoid the storage of N observation
# pairs and perform MVTMLE0 "in loops":
{
N <- n*(n-1)/2

Psi <- LocalPsi(Xs,nu,n,p,N)
Tmp <- eigen(Psi,symmetric=TRUE)
nG <- sqrt(sum((1 - Tmp\$values)^2))
iter <- 0
while (nG > delta && iter < 1000)
{
iter <- iter + 1
B <- B %*% Tmp\$vectors
Xs <- Xs %*% Tmp\$vectors
a <- LocalA(Xs,nu,n,p,N,Tmp\$values)

# Check whether the new matrix
# parameter is really better:
a.half <- a/2
Xs.new <- t(t(Xs)*exp(-a.half))
DL <- LocalDL(Xs,nu,n,p,N,Xs.new,a)
DL0 <- sum(a*(1 - Tmp\$values))/4
if (steps)
{
print(cbind(Before.iteration=iter,
Diff.M2LL.with.PN=DL))
}
if (DL <= DL0)
{
B <- t(t(B) * exp(a.half))
Xs <- Xs.new
} else
# perform a step of the classical
# fixed-point algorithm:
{
sqrtTmpValues <- sqrt(Tmp\$values)
B <- t(t(B) * sqrtTmpValues)
Xs <- t(t(Xs) / sqrtTmpValues)
}
Psi <- LocalPsi(Xs,nu,n,p,N)
Tmp <- eigen(Psi,symmetric=TRUE)
nG <- sqrt(sum((1 - Tmp\$values)^2))
}
if (steps)
{
print(cbind(After.iteration=iter,
}
}
S <- tcrossprod(B)
if (nu == 0)
{
S <- S / det(S)^(1/p)
}
rownames(S) <- colnames(S) <- Xnames
return(list(B=B,S=S,iter=iter))
}

LocalPsi <- function(X,nu,n,p,N)
{
Y <- X[n,] - X[n-1,]
Z <- Y/sqrt(nu + sum(Y^2))
Psi <- tcrossprod(Z)
for (i in 1:(n-2))
{
jj <- (i+1):n
Y <- t(t(X[jj,]) - X[i,])
denom <- nu + rowSums(Y^2)
Z <- Y/sqrt(denom)
Psi <- Psi + crossprod(Z)
}
Psi <- Psi * ((nu + p)/N)
return(Psi)
}

LocalA <- function(X,nu,n,p,N,evs)
{
Y <- X[n,] - X[n-1,]
Z <- Y^2/(nu + sum(Y^2))
Ht <- tcrossprod(Z)
for (i in 1:(n-2))
{
jj <- (i+1):n
Y <- t(t(X[jj,]) - X[i,])
denom <- nu + rowSums(Y^2)
Z <- Y^2/denom
Ht <- Ht + crossprod(Z)
}
Ht <- diag(evs) -
Ht * ((nu + p)/N) + (nu == 0)/p
a <- qr.solve(Ht,evs - 1)
return(a)
}

LocalDL <- function(X,nu,n,p,N,Xnew,a)
{
Y    <- X[n,]    - X[n-1,]
Ynew <- Xnew[n,] - Xnew[n-1,]
DL <- log((nu + sum(Ynew^2))/
(nu + sum(Y^2)))
for (i in 1:(n-2))
{
jj <- (i+1):n
Y    <- t(t(X[jj,])    - X[i,])
z    <- nu + rowSums(Y^2)
Ynew <- t(t(Xnew[jj,]) - Xnew[i,])
znew <- nu + rowSums(Ynew^2)
DL <- DL + sum(log(znew/z))
}
DL <- DL * ((nu + p)/N) + sum(a)
return(DL)
}
```

## Try the fastM package in your browser

Any scripts or data that you put into this service are public.

fastM documentation built on May 2, 2019, 4:01 a.m.