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

```# Auxiliary function:

SqrtS <- function(S)
# For a symmetric, positive definite matrix S,
# this procedure returns a matrix B such that
#    S = B %*% t(B).
{
res <- eigen(S,symmetric=TRUE)
B <- t(t(res\$vectors) * sqrt(res\$values))
return(B)
}

# Function to check how well a matrix parameter
# satisfies the fixed-point equation:

CheckS0 <- function(X,nu,S)
{
n <- dim(X)[1]
p <- dim(X)[2]
denom <- nu + rowSums(X*t(qr.solve(S,t(X))))
Z <- X/sqrt(denom)
S1 <- (nu + p) * crossprod(Z)/n
err1 <- norm(S-S1)
R <- SqrtS(S)
Y <- t(qr.solve(R,t(X)))
denom <- nu + rowSums(Y^2)
Z <- Y/sqrt(denom)
S2 <- (nu + p) * crossprod(Z)/n
err2 <- norm(S2 - diag(rep(1,p)))
return(c(err1=err1,err2=err2))
}

# Algorithm for the location-scatter problem:

MVTMLEr <- function(X,nu=1,delta=10^(-7))
{
if (nu < 1)
{
print("Need nu >= 1 degress of freedom!")
print("Use nu == 1 now!")
nu <- 1
}
Xnames <- colnames(X)
X <- as.matrix(X)
n <- dim(X)[1]
p <- dim(X)[2]
Y <- cbind(X,rep(1,n))
G <- MVTMLE0r(Y,nu-1,delta)\$S
G <- G / G[p+1,p+1]
mu.hat <- G[1:p,p+1]
names(mu.hat) <- Xnames
Sigma.hat <- G[1:p,1:p] - tcrossprod(mu.hat)
rownames(Sigma.hat) <-
colnames(Sigma.hat) <- Xnames
Shape.hat <- Sigma.hat / det(Sigma.hat)^(1/p)
return(list(mu.hat=mu.hat,
Sigma.hat=Sigma.hat,
Shape.hat=Shape.hat))
}
```

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.