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.
Best Answer
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") Similar Posts:
- Solved – Ellipse formula from points
- Solved – the proper way of calculating the kernel density estimate from geographical coordinates
- Solved – Condition number calculation in R
- Solved – the 95% confidence interval / ellipse in a PCA plot telling me
- Solved – How to get ellipse region from bivariate normal distributed data