[R] Cross-checking a custom function for separability indices

[Nikos Alexandris]

Hi list!

I have prepared a custom function (below) in order to calculate separability 
indices (Divergence, Bhattacharyya, Jeffries-Matusita, Transformed divergene) 
between two samples of (spectral land cover) classes.

I need help to cross-compare results to verify that it works as expected 
(since I don't know of any other foss-tool that will give me quickly some 
results).

Does anybody use another mean/tool/function (within or out of R) to calculate 
these distances?


For example,

# two samples
sample.1 <- c (1362, 1411, 1457, 1735, 1621, 1621, 1791, 1863, 1863, 1838)
sample.2 <- c (1354, 1458, 1458, 1458, 1550, 1145, 1428, 1573, 1573, 1657)

# running the custom function (below)
separability.measures ( sample.1 , sample.2 )

Divergence:             1.5658 
Bhattacharryya:         0.1683 
Jeffries-Matusita:      0.3098 
Transformed divergence: 0.3555

Thanks in advance, Nikos



# --- Code -------------------------------------------------------------------
# Custom function for various Separability Measures
  # by Nikos Alexandris, Freiburg, 8.04.2010


# ( based on Divergence and Jeffries-Matusita, requires input variables 
"as.matrices" )
separability.measures <- function ( Vector.1 , Vector.2 ) {


    # convert vectors to matrices in case they are not
    Matrix.1 <- as.matrix (Vector.1)
    Matrix.2 <- as.matrix (Vector.2)


    # define means
    mean.Matrix.1 <- mean ( Matrix.1 )
    mean.Matrix.2 <- mean ( Matrix.2 )


    # define difference of means
    mean.difference <- mean.Matrix.1 - mean.Matrix.2


    # define covariances for supplied matrices
    cv.Matrix.1 <- cov ( Matrix.1 )
    cv.Matrix.2 <- cov ( Matrix.2 )


    # define the halfsum of cv's as "p"
    p <- ( cv.Matrix.1 + cv.Matrix.2 ) / 2


# --%<------------------------------------------------------------------------
    # calculate the Bhattacharryya index
    bh.distance <- 0.125 *
                         t ( mean.difference ) *
                         p^ ( -1 ) *
                         mean.difference +
                         0.5 * log10 (
                                     det ( p ) / sqrt (
                                                      det ( cv.Matrix.1 ) *
                                                      det ( cv.Matrix.2 )
                                                      )
                                     )


# --%<------------------------------------------------------------------------
    # calculate Jeffries-Matusita
    jm.distance <- 2 * ( 1 - exp ( - bh.distance ) )


# --%<------------------------------------------------------------------------
    # calculate the divergence
      # trace (is the sum of the diagonal elements) of a square matrix
      trace.of.matrix <- function ( SquareMatrix ) {
                     sum ( diag ( SquareMatrix ) ) }

      # term 1
      divergence.term.1 <- 1/2 *
                 trace.of.matrix (
                                  ( cv.Matrix.1 - cv.Matrix.2 ) *
                                  ( cv.Matrix.2^ (-1) - cv.Matrix.1^ (-1) )
                                 )

      # term 2
      divergence.term.2 <- 1/2 *
                 trace.of.matrix (
                                  ( cv.Matrix.1^ (-1) + cv.Matrix.2^ (-1) ) *
                                  ( mean.Matrix.1 - mean.Matrix.2 ) *
                                  t ( mean.Matrix.1 - mean.Matrix.2 )
                                 )

      # divergence
      divergence <- divergence.term.1 + divergence.term.2


# --%<------------------------------------------------------------------------
    # and the transformed divergence
    transformed.divergence  <- 2 * ( 1 - exp ( - ( divergence / 8 ) ) )


# --%<------------------------------------------------------------------------
    # print results --- look at "separability.measures.pp"

    # get "column" names (hopefully they exist) --- to use for/with... ??? ---
    name.of.Matrix.1 <- colnames ( Matrix.1 )
    name.of.Matrix.2 <- colnames ( Matrix.2 )


    # create new objects to be used with... ?
    assign ( "divergence.vector"             , divergence             ,
 envir = .GlobalEnv)

    assign ( "jm.distance.vector"            , jm.distance            ,
 envir = .GlobalEnv)

    assign ( "bh.distance.vector"            , bh.distance            ,
 envir = .GlobalEnv )

    assign ( "transformed.divergence.vector" , transformed.divergence ,
 envir = .GlobalEnv )


    # print some message --- ?
    # a "just do it" solution using cat()
    cat ( paste ( "Divergence:             " ,
 round ( divergence , digits = 4 ) , sep = "" ) , "\n")

    cat ( paste ( "Bhattacharryya:         " ,
 round ( bh.distance , digits = 4 ) , sep = "" ) , "\n")

    cat ( paste ( "Jeffries-Matusita:      " ,
 round ( jm.distance , digits = 4 ) , sep = "" ) , "\n")

    cat ( paste ( "Transformed divergence: " ,

 round ( transformed.divergence , digits = 4 ) , sep = "" ) , "\n" )
    cat ("\n")


}
# --- Code -------------------------------------------------------------------


# --%<------------------------------------------------------------------------
# Sources

  # Richards, John A., Jia, X., "Remote Sensing Digital Image Analysis: ...", 
Springer
  # Wikipedia