R/robust.R

In rotations: Working with Rotation Data

#' Measure of Discord
#' 
#' This function computes a measure of discord for a sample of random rotations.
#' The larger the statistic value the less likely it is the corresponding observation
#' was generated by the same mechanism the rest of the data as generated by.
#' It can be used to test for outliers in SO(3) by comparing it to an F distribution
#' with 3,3(n-2) df for the Cayley or matrix Fisher distributions or to an F
#' distribution with 1,n-2 df for the von Mises Fisher distribution.
#' 
#' @param x The sample of random rotations
#' @param type To specify if "intrinsic" or "extrinsic" approach should be used to compute the statistic
#' @param t If test blocs then the bloc size, set to 1 by default
#' @param obs integer vector specifying which observation(s) to compute the measure of discord for
#' @return The Hi statistic for each group of size t is returned.  If \code{t>1} then which observations 
#' that define each group of size \code{t} is returned as well.
#' @export
#' @examples
#' #Compute the measures of discord for a sample from the Cayley distribution
#' # Intrinsic examples are commented out but are below if you're interested
#' 
#' Rss <- ruars(20,rcayley,kappa=1)
#' Hi <- discord(Rss, type='intrinsic')
#' He <- discord(Rss, type='extrinsic')
#' 
#' #Compare to the theoretical F distribution
#' OrdHi <- sort(Hi)
#' OrdHe <- sort(He)
#' 
#' par(mfrow=c(1,2))
#' plot(ecdf(OrdHi),main='Intrinsic',xlim=range(c(OrdHi,OrdHe)))
#' lines(OrdHi,pf(OrdHi,3,3*(length(OrdHi)-2)))
#' 
#' plot(ecdf(OrdHe),main='Extrinsic',xlim=range(c(OrdHi,OrdHe)))
#' lines(OrdHi,pf(OrdHi,3,3*(length(OrdHe)-2)))
#' layout(1)

discord<-function(x, type, t=1L, obs=1:nrow(x)){
  #Compute the statistic proposed by Best and Fisher (1986) that is a function of the largest eigenvalue
  #when observation i was removed
  #Written for quaternions, so if SO3 is given, make them quaternions
  
  Qs<-as.Q4(x)
  type <- try(match.arg(type,c("intrinsic", "extrinsic")),silent=T)
  if ("try-error"%in%class(type))
    stop("type needs to be one of 'intrinsic' or 'extrinsic'.")
  
  if(any(obs>nrow(x))||any(obs<0))
    stop("obs must be between 1 and nrow(x)")
  
  if(t>1){
    if(type=='intrinsic'){
      warning("intrinsic en bloc testing isn't available yet")
    }
    return(HnBlocCpp(Qs,t))
    
  }else if(t<1){
    
    stop("t must be an at least 1")
    
  }
  
  #If t==1 then an Hi statistic is computed for each observation
  if(type=='extrinsic'){
    
    Hn <- as.vector(HnCpp(Qs))[obs]
    
  }else if(type=='intrinsic'){
    
    Hn <- as.vector(HnCppIntrinsic(Qs))[obs]
    
  }else{
    
    stop("Please choose extrinsic or intrinsic type.")
    
  }
  
  if(any(Hn<0)){
    HnSub<-as.vector(HnCpp(Qs))
    ToReplace<-which(Hn<=0)
    Hn[ToReplace]<-HnSub[ToReplace]
    warning(paste(length(ToReplace)," Hn values were negative.  They were replaced by the extrinsic approximations."))
  }
  
  return(Hn)
}

HnBlocCpp<-function(Qs,t){
  #Compute the Hn statistic when each possible set of t observations is deleted
  Qs<-as.Q4(Qs)
  n<-nrow(Qs)
  groups <- utils::combn(n,t)
  
  Hnia <- HnCppBloc(Qs,groups)
  
  return(list(Groups=groups,Hi=Hnia))
  
}