# Solved – 95% confidence ellipse of Lenth’s maximum likelihood estimation

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(spi/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= -sumlimits_{i=1}^n(y-y_i)[s_i(x-x_i)-c_i(y-y_i)]/d^3_i = 0\$

\$L_y= sumlimits_{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.

Contents

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!
``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") ``