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(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.
 A: 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")

A: 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.
