I wish to perform Jeffries-Matusita distance on 14 spectral bands. Is there anyone who can help with how it is done in R?
-
2Questions asking for code are off-topic on CV (see our [help center](http://stats.stackexchange.com/help/on-topic)). They can be on-topic on [Stack Overflow](http://stackoverflow.com/), but this wouldn't be. You should take their [tour](http://stackoverflow.com/tour), & read their [guidelines for asking questions](http://stackoverflow.com/help/asking) there. Then you may be able to ask an updated version of this question there. – gung - Reinstate Monica Jul 09 '14 at 13:35
-
2This question appears to be off-topic because it is about asking for code. – gung - Reinstate Monica Jul 09 '14 at 13:35
-
I have edited the question.cheers. – tribalhouse49778 Jul 09 '14 at 13:48
-
1This is still asking for code, & is still off-topic here & on SO. – gung - Reinstate Monica Jul 09 '14 at 14:00
-
If you need general formulas or explanations, those would be on-topic here--and it is likely that any answers would illustrate them with some `R` code ;-). – whuber Jul 09 '14 at 14:01
-
The spelling Jeffries seems common in the literature, but my guess is that it should be Jeffreys. @whuber – Nick Cox Oct 14 '19 at 12:33
-
2@Nick Good call. I can trace the history to H. Jeffreys (1946), "An Invariant for the Prior Probability in Estimation Problems," cited by Wacker, A. G. and D. A. Landgrebe (1972), "Minimum Distance Classification in Remote Sensing." A year later NASA researchers misspelled his name as "Jeffries" in Swain, P.H. (1973), "Pattern recognition: a basis for remote sensing data analysis" and it looks like everyone thereafter used the misspelling. (Which is strong evidence nobody was reading the original--they just kept citing the citation.) – whuber Oct 14 '19 at 12:47
-
2Sir Harold himself told me (c. 1976) that he disliked the form Jeffreys' as it would all too easily morph into Jeffrey's but he said nothing about Jeffries, which is just wrong. Astronomers are especially prone to conflate Harold Jeffreys and William H. Jefferys. – Nick Cox Oct 14 '19 at 12:53
1 Answers
As Nick Cox points out in a comment to the question, the Jeffries-Matusita distance should be called the Jeffreys-Matusita distance due to its origin in the work of Harold Jeffreys.
Whatever you choose to call it, distance between two continuous distributions $F$ and $G$ with PDFs $f$ and $g$ is a re-expression of their Bhattacharyya distance:
$$D_{JM}(F,G) = \sqrt{2(1 - \exp(-D_B(F,G))}$$
where
$$D_B(F,G) = \int_\mathbb{R}\sqrt{f(x)g(x)}dx.$$
For nondegenerate multivariate Normal distributions with means $\mu_i$ and covariance matrices $\Sigma_i$ this works out to
$$\eqalign{ &D_B(\mathcal{N}(\mu_F,\Sigma_F), \mathcal{N}(\mu_G,\Sigma_G)) = \\ &\frac{1}{8}D_\Sigma(\mu_F,\mu_G) + \frac{1}{2}\left(\log(|\Sigma|) - \log(|\Sigma_F|)/2 - \log(|\Sigma_G|)/2\right) }$$
where $\Sigma=(\Sigma_F+\Sigma_G)/2$ is the component-wise average of the two covariance matrices, $|\cdot|$ denotes the determinant, and $D_\Sigma$ is the Mahalanobis distance with respect to $\Sigma$,
$$D_\Sigma(\mu_F, \mu_G) = (\mu_F - \mu_G)^\prime \Sigma^{-1}(\mu_F - \mu_G).$$
I presume that the distance is desired between two "signatures." These are estimated multivariate normal distributions of the vectors $(x_1,x_2,\ldots,x_{14})$ associated with distinct groups of pixels in the image. (The components of these vectors are the "bands".)
Implementation considerations--all platforms
There are some important implementation details to note, especially for those who wish to generalize this to larger multispectral data, which can have hundreds or even thousands of bands. In these cases the covariance matrices become huge and their determinants can easily overflow double precision floating point numbers. In fact, it's usually not a good idea to compute determinants of large matrices in the first place, but if you have to, an efficient accurate way is to reduce each matrix to a triangular form. Once that is done, the determinant is found as the product of the diagonal elements. Equivalently, the log determinant (which is involved in $D_{JM}$) can be found by summing the logs of the diagonal elements. This avoids over- and underflow.
It also is not necessary to invert $\Sigma$, since all that is needed is the product $\nu=\Sigma^{-1}(\mu_F-\mu_G)$. Find this instead by solving the system of equations
$$\Sigma(\nu)= \mu_F-\mu_G.$$
That will be much faster and a little more precise.
It can be flexible and convenient to implement all three distances mentioned here, because this assists testing and provides additional functionality. Thus, using R
to illustrate, we might code something like these functions:
#
# Compute the Mahalanobis distance between two vectors.
#
mahalanobis <- function(m1, m2, sigma) {m <- m1 - m2; m %*% solve(sigma, m)}
#
# Compute the Bhattacharyya distance between two multivariate normal distributions
# given by their means and covariance matrices.
#
bhattacharyya <- function(m1, s1, m2, s2) {
d <- function(u) determinant(u, logarithm=TRUE)$modulus # Log determinant of matrix u
s <- (s1 + s2)/2 # mean covariance matrix
mahalanobis(m1, m2, s)/8 + (d(s) - d(s1)/2 - d(s2)/2)/2
}
#
# Re-express the Bhattacharyya distance as the Jeffries-Matusita distance.
#
jeffries.matusita <- function(...) sqrt(2*(1-exp(-bhattacharyya(...))))
That is simple, straightforward, fast, and readily portable to most computing platforms.
Illustration in R
An example of the use of these functions in R
follows. It generates sample multi-spectral data for a given set of classes, representing them as an array in which each band's values appear in the columns together with a parallel array identifying the class of each row (aka "pixel") with an index $1, 2, \ldots$. It computes the class "signatures" (the estimated multivariate normal parameters). The signature data structure is a list consisting of the mean vector and covariance matrix; the class signatures are collected into one large list. It then applies $D_{JM}$ to all pairs of distinct distance signatures in that list to produce a distance matrix indexed by the classes.
#
# Generate sets of data for d bands aggregated into classes.
#
d <- 14 # Number of bands
n <- rep(1000, 5) # Class pixel counts
n.class <- length(n) # Number of classes
require(MASS) # For generating multivariate normals
set.seed(17) # Allows reproducible results
# Create random mean vectors for the classes
mu <- round(matrix(rnorm(d*n.class, 128, 1), ncol=n.class, byrow=TRUE), 0)
# Initialize the data structure {x, classes}
x <- matrix(double(), ncol=d, nrow=0)
classes <- integer()
# Generate the random data
for (i in 1:n.class) {
# Create a random valid covariance matrix for this class
f <- svd(matrix(rnorm(d^2), ncol=d))
sigma <- t(f$v) %*% diag(rep(10, d)) %*% f$v
# Generate d-variate normals
x <- rbind(x, mvrnorm(n[i], mu[, i], sigma))
classes <- c(classes, rep(i, n[i]))
}
#------------------------------------------------------------------------------------#
#
# Given a set of bands (as the columns of x) and a grouping variable in `class`,
# compute the class means and covariance matrices (the "signatures").
#
classes.stats <- by(x, classes, function(y) list(mean=apply(y, 2, mean), cov=cov(y)))
#------------------------------------------------------------------------------------#
#
# Compute the J-M distances between the classes.
#
distances <- matrix(0.0, n.class, n.class)
for (i in 2:n.class) {
m1 <- classes.stats[[i]]$mean; s1 <- classes.stats[[i]]$cov
for (j in 1:(i-1)) {
m2 <- classes.stats[[j]]$mean; s2 <- classes.stats[[j]]$cov
distances[i,j] <- distances[j,i] <- jeffries.matusita(m1,s1,m2,s2)
}
}
print(distances)
-
I am encountering an error on the last function. Its saying "Error: could not find function "jeffries.matusita"" – tribalhouse49778 Jul 10 '14 at 12:03
-
It appears you did not implement the first block of code in which that function is defined. – whuber Jul 10 '14 at 13:06