R/marpca.sub.R

In musto101/wilcox_R: What the package does (short line)

marpca.sub <-
function(x,p=ncol(x)-1,N1=3,N2=2,tol=.001,B=NULL,a=NULL,
LSCALE=TRUE){
#
# Marrona (2005, Technometrics, 47, 264-273) robust PCA
#
# Note: setting
# p=0 causes B to be the identity matrix, which is used in the case
# p=ncol(x) to estimate proportion of unexplained variance.
#
wt.cov<-NULL
if(!is.null(B)){
B<-as.matrix(B)
if(ncol(B)==1)B<-t(B)
}
n<-nrow(x)
m<-ncol(x)
q<-m-p
if(q<0)stop("p and q should have values between 1 and ncol(x)")
hval<-floor((n + m - q + 2)/2)
DEL<-Inf
sig0<-Inf
if(is.null(B)){
if(p>0 && p<m){
B<-matrix(runif(q*m),nrow=q,ncol=m)
B<-t(ortho(t(B)))
}
if(p==0)B<-diag(rep(1,m))
}
it<-1
Bx<-matrix(NA,ncol=q,nrow=n) #q by n
for(i in 1:n)Bx[i,]<-B%*%as.matrix(x[i,]) # q by 1
#Bx<-as.matrix(elimna(Bx))
if(is.null(a))a<-apply(Bx,2,FUN="median") # Initial a (coordinatewise medians)
while(it<N1+N2 && DEL>tol){
r<-NA
for(i in 1:n)r[i]<-sum(Bx[i,]-a)^2
if(LSCALE)sig<-lscale(r,m,q)
if(!LSCALE){
delta<-delta<-(n-m+q-1)/(2*n)
sig<-mscale(r,delta)
}
DEL<-1-sig/sig0
sig0<-sig
ord.r<-order(r)
w<-rep(0,n)
w[ord.r[1:hval]]<-1
xx<-x
for(i in 1:n)xx[i,]<-x[i,]*w[i]
mu<-apply(xx,2,FUN="sum")/sum(w) #m by 1 locations
Cmat<-matrix(0,nrow=m,ncol=m)
for(i in 1:n){
temp<-w[i]*as.matrix(x[i,]-mu)%*%t(as.matrix(x[i,]-mu))
Cmat<-Cmat+temp
}
wt.cov<-Cmat/sum(w)
if(it>N1){
temp<-eigen(wt.cov)
ord.eig<-order(temp$values)
for(iq in 1:q)B[iq,]<-temp$vectors[,ord.eig[iq]]
}
a<-B%*%mu
it<-it+1
}
list(B=B,a=a,var.op=sig,mu=mu,wt.cov=wt.cov)
}