95% confidence ellipse of Lenth's maximum likelihood estimation

Question

I'm attempting to write a simple mobile app to help colleagues triangulate the signal of transmitter-fitted animals in the field. I've been using Russell Lenth's 1981 method of Maximum Likelihood Estimation (PDF: On Finding the Source of a Signal)

Basically, I start with 3+ bearings and the locations (in UTMs) they were taken from. I can iterate through solving equality 2.6 in the PDF to get the estimated location easily enough, but now I'd like to be able to calculate the area of a 95% confidence ellipse of the estimated location ($x$, $y$). Ideally, I would also like to be able to draw this confidence ellipse on a map.

However, I'm not a statistics person and despite looking up the explanations for covariance estimates and concentration parameters, I'm a bit stuck on how to tackle this problem.

From sifting through a few books on wildlife tracking, I suspect that I need to know the standard deviation of the bearings (s) (figured out through system testing, I assume?) and then I can find $\kappa$ using the following:

$A = exp[-1/2(s\pi/180)^2]$

$\kappa^{-1} = 2(1-A) + [(1-A)^{2}(0.48794-0.82905A-1.3915A^{2})]/A $

Am I on the right track so far? And if so, any hints on what my next step is? I'm guessing it has something to do with $Q$ in the PDF? Any help is greatly appreciated.

---

Example: Say I start with the following:

Point 1: 554045 Easting, 2813867 Northing, Bearing: 300$^{\circ}$

Point 2: 553355 Easting, 2813873 Northing, Bearing: 20$^{\circ}$

Point 3: 553207 Easting, 2814120 Northing, Bearing: 90$^{\circ}$

I believe the likelihood described in the paper is:

$L_x= -\sum\limits_{i=1}^n(y-y_i)[s_i(x-x_i)-c_i(y-y_i)]/d^3_i = 0$

$L_y= \sum\limits_{i=1}^n(x-x_i)[s_i(x-x_i)-c_i(y-y_i)]/d^3_i = 0$

With

$s_i = \sin \theta_i$ ($\theta$ is the bearing in radian measure and in the mathematical sense)

$c_i = \cos \theta_i$

$d_i = [(x-x_i)^{2}+(y-y_i)^{2}]^{1/2}$

I solved Equality 2.6 in the PDF (sorry, my LaTeX isn't that strong) by first ignoring the asterisks to find a starting estimate, and then using that estimate to find interim estimates of $d$, $s^*_i$, and $c^*_i$ (defined below), which were plugged into Equation 2.6 and then this was repeated until the ($x$, $y$) findings converged to a single UTM point.

$s^*_i = (y-y_i)/d^3_i$

$c^*_i = (x-x_i)/d^3_i$

$d_i = [(x-x_i)^{2}+(y-y_i)^{2}]^{1/2}$

$z_i = s_ix_i - c_iy_i$

Using this method, the bearings in the Example would converge at 553455 Northing, 2814131 Easting. I would like to now figure out what the area of the 95% confidence ellipse of the estimated location is, as well as how it would be drawn on a map.

---

Edit: Whoops, apparently I haven't been getting email notifications on this post. Apologies for any ignored comments. I took my problem to a colleague that pointed me to this PDF which explains error ellipses in a much more digestible way than I've seen before. I've only run through a few examples so far, but as far as I can tell, the ellipses look the way I would expect them to.

Answer 1

Here is an R implementation. With it I find that the algorithm does not always converge; sometimes, even when it does, it produces negative covariances; and it tends to be optimistic concerning the coverage of the confidence ellipse (about half the 95% confidence ellipses do not cover the true values). That could be due to bugs on my part, so check and test carefully!

In this simulation result the true location is in blue, the estimated location based on 16 very noisy bearings is in red, the bearings from the reference points are gray rays, and the 95% confidence ellipse is shown centered around the estimate.

library(MASS) # ginv()
lenth <- function(xy, theta, x.hat=NULL, i.max=1000, threshold=10.0^(-6)) {
  norm2 <- function(v) sum(v*v)
  #
  # Pseudo-MLE.
  #
  improve <- function(x, cs, xy){
    #
    # Improves an initial estimate `x` based on bearings as represented
    # by a 2 x * matrix of cosines and sines `cs` relative to locations `xy` 
    # (as a 2 by * matrix).
    #
    if (is.null(x)) {
      cs.star <- cs                                          # Starting estimate
    }
    else {
      xy.relative <- apply(xy, 2, function(v) x-v)           # vectors towards x
      d <- apply(xy.relative, 2, function(v) norm2(v)^(3/2)) # cubed distances
      cs.star <- xy.relative / d                             # direction cosines
    }
    z <- cs[2,]*xy[1,] - cs[1,]*xy[2,]                     
    a <- cs.star %*% t(cs) * matrix(c(-1,-1,1,1), nrow=2)    # Eq 2.6, lhs
    b <- cs.star %*% z                                       # Eq 2.6, rhs
    u <- solve(a,b)                                          # Solution
    c(u[2], u[1])
  }
  #
  # Convert bearings into direction cosines and sines.
  #
  cs <- sapply(theta, function(x) c(cos(x), sin(x)))
  #
  # Iterate until convergence.
  #
  if (is.null(x.hat)) x.hat <- improve(x.hat, cs, xy)
  repeat {
    x.new <- improve(x.hat, cs, xy)
    eps <- norm2(x.hat-x.new) / sqrt(max(apply(xy, 2, norm2)))
    lines(t(cbind(x.hat, x.new)), pch=2, col="Red")
    x.hat <- x.new
    i.max <- i.max - 1
    if (eps <= threshold) break
    if (i.max < 0) stop("Convergence failed.")
  }
  #
  # Estimate kappa.
  #
  mu.hat <- apply(x.hat - xy, 2, function(u) atan2(u[2], u[1]))
  c.bar <- mean(cos(theta - mu.hat))
  kappa.inv <- 2*(1-c.bar) + (1-c.bar)^2 * (0.48794/c.bar - 0.82905 - 1.3915*c.bar)
  kappa <- 1/kappa.inv
  #
  # Estimate the covariance matrix for x.hat.
  #
  xy.relative <- apply(xy, 2, function(v) x.hat - v)
  d <- apply(xy.relative, 2, function(v) norm2(v)^(3/2)) 
  cs.star <- xy.relative / d  
  q.hat <- cs.star %*% t(cs)
  q.hat <- (q.hat + t(q.hat))/2
  q.hat <- q.hat * matrix(c(1,-1,-1,1), nrow=2)
  q.hat <- q.hat[2:1, 2:1]
  q.hat <- ginv(q.hat, tol=10^(-16)) * kappa.inv
  x.se <- sqrt(q.hat[1,1])
  y.se <- sqrt(q.hat[2,2])

  list(x.hat=x.hat, cov=q.hat, se=c(x.se, y.se), kappa=kappa)
}
#--------------------------------------------------------------------------------------#
set.seed(17)
#
# Target point.
#
x.0 <- c(1, 1)
#
# Observation points.
#
n <- 16
xy <- matrix(rnorm(2*n), nrow=2)
cs <- sapply(theta, function(x) c(cos(x), sin(x)))
#
# Simulate the observed bearings (in radians)
#
xy.fuzzy <- matrix(rnorm(2*n, sd=0.2), nrow=2) + x.0
theta <- apply(xy.fuzzy - xy, 2, function(u) atan2(u[2], u[1]))
#
# Obtain the MLE.
#
fit <- lenth(xy, theta)
#
# Plot all points.
#
cs <- sapply(theta, function(x) c(cos(x), sin(x)))
par(mfrow=c(1,1))
plot(t(cbind(xy, fit$x.hat, x.0)), asp=1, type="n", xlab="x", ylab="y")
temp <- sapply(1:length(theta), 
           function(i) lines(t(cbind(xy[,i], xy[,i]+9*cs[,i])), col="Gray"))
points(t(xy))
points(x.0[1], x.0[2], pch=19, col="Blue") 
points(t(fit$x.hat), pch=19, col="Red")
#
# Plot the 95% confidence ellipse.
#
lines(ellipse(fit$cov, centre=fit$x.hat, level=0.95), col="#AA4444")

Answer 2

The link to Lenth's paper no longer works.

I have been working on a very similar application in Java. I recently translated the Fortran code for TRIANG (which utilizes Lenth's LME, and the source code is in the public domain) into Java. You may have an error in your code for the position estimate, although your interpretations of Lenth's equations look correct (but I did not go through each thoroughly). Using your example in my code, I get for a position estimate: Xcoor: 553650.1187572011 Ycoor: 2813504.348551782 iterations: 7 I get the same result with the TRIANG program and the LOCATE III program. To solve the covariance matrix I use the Blas_j.java, and linpack.java classes, also available on the internet. assuming a standard deviation of 2.5, I get an error ellipse: area of error ellipse: 2402.3515125037443 and a covariance matrix of: cv1,1: 168.93858246339704 cs1,2: 31.635641721682703 cv2,1: 31.635641721682703 cv2,2: 102.39381571783277

Now, my problem is similar to yours. I am trying to solve the 95% error ellipse. Given a prior SD estimate, I can confirm a correct ellipse area, but when I try to estimate the SD from my data, I get a different answer from the two programs listed above.

To plot the ellipse, I believe you need the eigenvalue and eigenvector from the above matrix to plot the ellipse.

EDIT: Oops. You are correct! After re-inputting your data I get the same values that you do - (Xcoor: 553455.1374349217 Ycoor: 2814130.724793033). And a variance matrix of cv1,1: 383.17463646469497; cv1,2: 71.75374233212538; cv2,1: 71.75374233212538; cv2,2: 232.24246671073084. So I think your position equations are working well. The errorEllipses PDF is very useful and practical guide to error ellipses. You might also be interested in this (warning, large download): http://demonstrations.wolfram.com/EllipseRepresentingTheConfidenceRegionOfACovarianceMatrix/ , which helped me to visualize how the CV matrix changes the shape and orientation of the error ellipse. However, even with all of this information, my error ellipses come out seemingly too small.

Answer 3

Thanks for taking this issue up. I was about to start trying to code this myself when I came across this thread!

I know the focus here is on calculating 95% CI ellipses, but the fact that the xy location you are estimating is not within the ellipse is not unusual.

Chp 3 in this book (http://books.google.com/books/about/Radio_Tracking_and_Animal_Populations.html?id=PcMd0fl55i8C) discusses how error polygons, ellipses, etc. frequently do not include the "real" location, as they are estimates of precision, not accuracy. They report on a beacon study they did; confidence ellipses only contained the "real" point about 60-65% of the time. I know it's not helpful, but it seems like this issue is not critical to wildlife studies at least.

Drawing a circle with, radius = linear error estimates, appears to be more accurate, but still not at that 95% level

Also, check out White and Garrot, 1990 "Analysis of Wildlife Radio-Tracking Data". They give SAS Macro coding for the MLE estimate and error ellipse (http://www.amazon.com/Analysis-Wildlife-Radio-Tracking-Data-White/dp/0127467254)

I was able to get both texts as pdf from my university library. Let me know if you have trouble finding them. Sorry if this is stuff you already knew.