tbs <-
function(x,eps=1e-3,maxiter=20,r=.45,alpha=.05,init.est=cov.mcd){
# Rocke's contrained s-estimator
#
# r=.45 is the breakdown point
# alpha=.05 is the asymptotic rejection probability.
#
library(MASS)
#if(!is.matrix(x))stop("x should be a matrix with two or more columns")
x<-elimna(x)
x=as.matrix(x)
n <- nrow(x)
p <- ncol(x)
LIST=F
if(p==1){
LIST=T
p=2
x=cbind(x,rnorm(nrow(x)))
# Yes, this code is odd, but for moment easiest way of handling p=1
}
temp<-init.est(x)
# very poor outside rate per obs under normality.
t1<-temp$center
s<-temp$cov
#if(p==1)stop("x should be a matrix with two or more columns")
c1M<-cgen.bt(n,p,r,alpha,asymp=FALSE)
c1<-c1M$c1
if(c1==0)c1<-.001 #Otherwise get division by zero
M<-c1M$M
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,center=t1,cov=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)
mnorm <- sqrt(as.vector(t.old) %*% as.vector(t.old))
snorm <- eigen(s.old)$values[1]
crit1 <- max(abs(t1 - t.old))
# crit <- max(crit1,crit2)
crit <- max(abs(w1d-wt.old))/max(w1d)
iter <- iter+1
}
if(LIST){
v1=t1[1]
v2=s[1,1]
return(list(center=v1,var=v2))
}
if(!LIST)return(list(center=t1,cov=s))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.