think the index is working now

ggobi · Mar 7, 2024 · faa90ce · faa90ce
1 parent 335541c
commit faa90ce
Show file tree

Hide file tree

Showing 4 changed files with 18 additions and 6 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -5,6 +5,8 @@
 * TeachingDemos removed as a Suggests, and replaced with aplpack for drawing Chernoff faces
 * addition of a pre-specified ellipse can be added to the 2D display
 * palette can now be a vector of values
+* a new projection pursuit index for finding anomalies relative to 
+a null variance-covariance matrix.
 
 # tourr 1.1.0 
 

diff --git a/R/anomaly-pursuit.r b/R/anomaly-pursuit.r
@@ -9,7 +9,7 @@ anomaly_index <- function() {
  function(mat, ell2d) {
 
  mat_tab <- #mean(mahal_dist(mat, ell2d))
- mean(sqrt(mahalanobis(mat, center=c(0,0), cov=ell2d)))
+ mean(mahalanobis(mat, center=c(0,0), cov=ell2d))
  }
 }
 
diff --git a/R/display-xy.r b/R/display-xy.r
@@ -153,11 +153,15 @@ display_xy <- function(center = TRUE, axes = "center", half_range = NULL,
  if (nrow(ellipse) == nrow(proj)) {
 
  # Project ellipse into 2D
+ # Notation in paper: ellipse=A, ellinv=A^(-1),
+ # e2=P^TA^(-1)P, ell2d=B
  evc <- eigen(ellipse*ellsize)
  ellinv <- (evc$vectors) %*% diag(evc$values) %*% t(evc$vectors)
  e2 <- t(proj) %*% ellinv %*% proj
  evc2 <- eigen(e2)
- ell2d <- (evc2$vectors) %*% diag(sqrt(evc2$values)) %*% t(evc2$vectors)
+ ell2d <- (evc2$vectors) %*% sqrt(diag(evc2$values)) %*% t(evc2$vectors)
+ e3 <- eigen(ell2d)
+ ell2dinv <- (e3$vectors) %*% diag(e3$values) %*% t(e3$vectors)
 
  # Compute the points on an ellipse
  sph <- geozoo::sphere.hollow(2, 200)$points
@@ -168,12 +172,14 @@ display_xy <- function(center = TRUE, axes = "center", half_range = NULL,
  sph2 <- sph2[order(sph2[,1], decreasing=T),]
  sph <- rbind(sph1, sph2)
  sph <- rbind(sph, sph[1,])
- sph2d <- sph%*%ell2d/half_range
+ sph2d <- sph%*%ell2dinv/half_range
 
  lines(sph2d)
 
  # Colour points outside the pD ellipse
- mdst <- sqrt(mahalanobis(data, center=rep(0, ncol(data)), cov=ellipse))
+ mdst <- mahalanobis(data,
+ center=rep(0, ncol(data)),
+ cov=ellipse)
  #mdst <- mahal_dist(data, ellipse)
  anomalies <- which(mdst > ellsize)
  points(x[anomalies,],

diff --git a/R/tour-guided-anomaly.r b/R/tour-guided-anomaly.r
@@ -46,7 +46,9 @@ guided_anomaly_tour <- function(index_f, d = 2, alpha = 0.5, cooling = 0.99,
 
  index <- function(proj) {
  # Check which observations are outside pD ellipse
- mdst <- sqrt(mahalanobis(data, center=rep(0, ncol(data)), cov=ellipse))
+ mdst <- mahalanobis(data,
+ center=rep(0, ncol(data)),
+ cov=ellipse)
  #mdst <- mahal_dist(data, ellipse)
  anomalies <- which(mdst > ellsize)
  stopifnot(length(anomalies) > 0)
@@ -56,7 +58,9 @@ guided_anomaly_tour <- function(index_f, d = 2, alpha = 0.5, cooling = 0.99,
  e2 <- t(proj) %*% ellinv %*% proj
  evc2 <- eigen(e2)
  ell2d <- (evc2$vectors) %*% diag(sqrt(evc2$values)) %*% t(evc2$vectors)
- index_f(as.matrix(data[anomalies,]) %*% proj, ell2d)
+ e3 <- eigen(ell2d)
+ ell2dinv <- (e3$vectors) %*% diag(e3$values) %*% t(e3$vectors)
+ index_f(as.matrix(data[anomalies,]) %*% proj, ell2dinv)
  }
 
  cur_index <- index(current)