R/GmedianCov.R

In Gmedian: Geometric Median, k-Medians Clustering and Robust Median PCA

GmedianCov <- function(X, init=NULL, nn=TRUE, scores=2, gamma=2, gc=2, alpha=0.75, nstart=1){
  ### Computation of the Geometric covariation matrix 
  ### with averaged stochastic gradient algorithms
  ### input : X   (n x p matrix, n observations, in dimension p)
  ### output : (geometric) median (1 x p numeric vector) and (geometric) median covariation matrix (p x p)  
  ### nn : non negative modification of the algorithm to get a non negaztive matrix estimation (TRUE by default)
  Gmed.est <- Gmedian(X,init=init,gamma=gamma,alpha=alpha,nstart=nstart)
  if (nn==FALSE) GMCM.est <- MedianCovMatRow_rcpp(X,Gmedian=Gmed.est,gamma=gc,alpha=alpha,nstart=nstart)
  else GMCM.est <- MedianCovMatRowP_rcpp(X,Gmedian=Gmed.est,gamma=gc,alpha=alpha,nstart=nstart)
  if (scores==0){ 
    return(list(median = Gmed.est,covmedian=GMCM.est))
  }
  else {
    ### Computation of the eigenvectors and scores 
    vectors <- RSpectra::eigs_sym(GMCM.est, scores)$vectors
    scores <- sweep(X,2,Gmed.est)%*%vectors
    sdev.rob <- apply(scores,2,scaleTau2) 
    return(list(median=Gmed.est,covmedian=GMCM.est,scores=scores,vectors=vectors,sdev=sdev.rob))
  }
}

WeiszfeldCov <- function(X, weights=NULL, scores=2, epsilon=1e-08, nitermax = 100){
  ### Computation of the Geometric covariation matrix 
  ### output : (geometric) median (1 x p numeric vector) and (geometric) median covariation matrix (p x p)  
  ### require library(RSpectra)
  X <- as.matrix(X)
  n <- nrow(X)
  if (is.null(weights)) poids <- rep(1,n) 
  else poids <- weights 
  Wmed.est <- Weiszfeld_rcpp(X,poids,epsilon=epsilon,nitermax=nitermax)
  WMCM.est <- MedianCovMatW_rcpp(X,Wmed.est$median,poids,epsilon=epsilon,nitermax=nitermax)
  if (scores==0){ 
    return(list(median = Wmed.est$median, covmedian=WMCM.est$median, iterm = Wmed.est$iter, itercov = WMCM.est$iter))
  }
  else {
    ### Computation of the eigenvectors and scores 
    vectors <- RSpectra::eigs_sym(WMCM.est$median, scores)$vectors
    vscores = sweep(X,2,Wmed.est$median)%*%vectors
    sdev.rob <- apply(vscores,2,scaleTau2) 
    return(list(median=Wmed.est$median, covmedian=WMCM.est$median, scores=vscores, sdev = sdev.rob, vectors=vectors, iterm = Wmed.est$iter, itercov = WMCM.est$iter))
  }
}