# What is the Maximum Likelihood Estimator for my nonlinear problem?

In the region $$0\leq x\leq 1$$ and $$0\leq y\leq 1$$ I have 1 radio broadcaster placed in $$(x_b,y_b)$$ and N receivers. The position of the broadcaster is not known. The position of each receiver however is known: $$(x_i,y_i)$$ as is the distance from the receiver to the broadcaster: $$d_i$$. There is some noise in the distance measurements that I assume to be normal distributed with mean=0 and the same variance for all receivers:

$$\sqrt{(x_i - x_b)^2+(y_i - y_b)^2} = d_i + \epsilon_i$$ where $$\epsilon_i \sim \mathcal{N}(0,\,\sigma^{2})$$

My goal is to estimate the most likely position of the radio broadcaster: $$(x_b,y_b)$$.

I want to do this by formulating this a minimization problem. So far I am minimizing: $$SE = \sum_i \left[\sqrt{(x_i - x_b)^2+(y_i - y_b)^2} - d_i\right]^2$$

I have read that for linear regression the Least Squares estimator is also the Maximum Likelihood Estimator under the assumption of normally distributed errors.

However is this also true for my nonlinear problem?

• The broadcaster-receiver distance, i.e. the value inside [...], can become negative, which does not fit your problem description. To stop this happening you would need non-normal errors. Commented Dec 15, 2023 at 8:56
• If you are interested in the maximum likelihood, it would be a good idea to state the (log-)likelihood Commented Dec 15, 2023 at 9:38
• A little tangent, but consider changing the distributional assumption of the $d_i$ from Normal to Rayleigh for a more realistic scenario Commented Dec 15, 2023 at 10:34
• I wonder if you're trying to reinvent the wheel. The wildlife literature has dealt with this for locating animals with radio collars for decades. See link.springer.com/article/10.1007/s10344-006-0076-9 and tandfonline.com/doi/abs/10.1080/00401706.1981.10486257.
– JimB
Commented Dec 15, 2023 at 17:15
• See stats.stackexchange.com/questions/625315/… for a similar problem where only the bearings to the source are known but not the distances; gis.stackexchange.com/a/2895/664 for an account of this kind of problem generally; and gis.stackexchange.com/a/40678/664 for the same problem with a relative error model (which might be more physically realistic in many applications).
– whuber
Commented Dec 15, 2023 at 19:25

I have read that for linear regression the Least Squares estimator is also the Maximum Likelihood Estimator under the assumption of normally distributed errors.

However is this also true for my nonlinear problem?

You are using

$$\sqrt{(x_i - x_b)^2+(y_i - y_b)^2} = d_i +\epsilon_i$$

but it is possibly better to state it as the following instead

$$\begin{array}{rcccl} d_i &=& \sqrt{(x_i - x_b)^2+(y_i - y_b)^2} &+& \epsilon_i \\ &=& \tilde{d}_i& + &\epsilon_i \end{array}$$

Then an observed distance $$d_{i}$$ is distributed according to the modeled mean distance $$\tilde{d}_i$$ plus some noise $$\epsilon_i$$.

Then the probability density (and the related likelihood) for the observations $$d_i$$ can be written as

$$\prod_{\forall i} \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( -\frac{(d_i-\tilde{d}_i)^2}{2\sigma^2} \right)$$

where the $$\tilde{d}_i = \sqrt{(x_i - x_b)^2+(y_i - y_b)^2}$$ can be modeled as functions of your unknown parameters $$y_b$$ and $$x_b$$.

Maximizing this likelihood is equivalent to minimizing the sum of squares $$\sum_{\forall i} (d_i-\tilde{d}_i)^2$$. It doesn't matter in what way $$\tilde{d}_i$$ is a function of the parameters.

The above is just to answer your direct question. Whether your model is sound (e.g. allowing negative distances) remains an open question.

Adding to Sextus Empiricus answer, here is a visualization of the likelihood in R, using the maximum likelihood solution for $$\sigma$$ from Ben's answer:

# plot circles around the respective points with the distances as radii
draw_circle <- function(x, y, r, col = "black", lty = 1, lwd = 1) {
theta <- seq(0, 2*pi, length.out = 100)
lines(x + r*cos(theta), y + r*sin(theta), col = col, lty = lty, lwd = lwd)
}

# define the likelihood function for p = c(x, y)
likelihood <- function(x, y, d_vec, p_vec) {
# true distances
d_true = sqrt((p_vec[,1] - x)^2 + (p_vec[,2] - y)^2)
# distance difference
d_diff = d_vec - d_true
# log-likelihood
log_likelihood = -sum(d_diff^2)
return(log_likelihood)
}

# plot the likelihood function for different values of xb and yb
plot_likelihood <- function(d_vec, p_vec, real_sigma, xlim, ylim,
n=101, pb) {
# define the grid
x_grid = seq(xlim[1], xlim[2], length.out = n)
y_grid = seq(ylim[1], ylim[2], length.out = n)
grid = expand.grid(x_grid, y_grid)
# compute the likelihood for each grid point
log_likelihood_grid = apply(grid, 1, function(x)
likelihood(x[1], x[2], d_vec, p_vec))
log_likelihood_grid = matrix(log_likelihood_grid, nrow = n,
ncol = n)
# add a point at the maximum
max_idx = which(log_likelihood_grid ==
max(log_likelihood_grid), arr.ind = TRUE)
xml = x_grid[max_idx[1]]
yml = y_grid[max_idx[2]]
# compute distances to the maximum
dml = sqrt((p_vec[,1] - xml)^2 + (p_vec[,2] - yml)^2)
# plug in the sigma^2 ml estimate
# I choose a biased lower variance correction, for this
# reason, this will only work with n>2 sensors
sigma2_ml = sum((d_vec - dml)^2)/(length(d_vec)-2)
log_likelihood_grid = -log(2*pi*sigma2_ml)/2 +
log_likelihood_grid/(2*sigma2_ml)
likelihood_grid =
exp(log_likelihood_grid)/sum(exp(log_likelihood_grid))
# set aspect ratio to 1
par(pty = "s")
# plot the likelihood
image(x_grid, y_grid, likelihood_grid, col = grey(seq(0, 1,
length.out = 100)), useRaster = TRUE, xlab = "x",
ylab = "y", xlim = c(0,1), ylim = c(0,1))
# plot p_vec
points(p_vec[,1], p_vec[,2], pch = 3, col = "cyan",
bg="black")
# draw circles around the points with the distances as radii
pd_vec = cbind(p_vec, d_vec)
apply(pd_vec, 1, function(x) draw_circle(x[1], x[2], x[3], "cyan"))
legend("topright", legend = c("Emitter", "Sensor", "ML",
"CI (99%)"), pch = c(16, 3, 16, NA),
col = c("red", "cyan", "blue", "blue"), bg = "black",
text.col = "white", box.col = "white",
lty=c(NA, NA, NA, 2), lwd=c(NA, NA, NA, 1))
# legend for the image
legend("bottomright", legend = c("Low likelihood",
"High likelihood"), fill = grey(
seq(0, 1, length.out = 100))[c(1, 100)],
bg = "black", text.col = "white",
box.col = "white", border="white")
ord_likeli = order(c(likelihood_grid), decreasing = TRUE)
reord_likeli= order(ord_likeli, decreasing = FALSE)
cumsum_prob = cumsum(likelihood_grid[ord_likeli])
cumsum_prob = cumsum_prob[reord_likeli]
cumsum_prob = t(matrix(cumsum_prob, nrow = n, ncol = n,
byrow = TRUE))
contour(x_grid, y_grid, cumsum_prob, add = TRUE, col = "blue",
levels=0.99, lty=2)
# plot the real position of the emitter
points(pb[1], pb[2], pch = 16, col = "red")
points(xml, yml, pch = 16, col = "blue")
# add title with estimate of sigma and the real sigma
title(bquote(paste(hat(sigma), " = ", .(round(sqrt(sigma2_ml),
3)), ", ", sigma, " = ", .(round(real_sigma, 3)))),
line = 0.5)
}


With randomly placed sensors and a large $$\sigma$$, to illustrate the problem better. Notice this is an approximate inference: the finer the grid is, the better it becomes.

n=5
p_vec = matrix(runif(2*n, 0, 1), ncol = 2)
pb = c(0.5, 0.5)
# real emmiter position
pb = c(0.5, sqrt(3)/6)

# distances, assumed to be normally distributed with equal
# variance around the true distance
d_true = sqrt((p_vec[,1] - pb[1])^2 + (p_vec[,2] - pb[2])^2)
sigma = 0.1
d_vec = d_true + rnorm(n, 0, sigma)

plot_likelihood(d_vec, p_vec, sigma, xlim = c(0, 1),
ylim = c(0, 1), pb=pb, n=1001)


• Those graphs might look nicer when you make the horizontal and vertical axes the same scale. Commented Dec 16, 2023 at 16:47
• See gis.stackexchange.com/a/40678/664 for an example of this kind of illustration. I would draw your attention to (a) the 1:1 aspect ratio, as pointed out by @Sextus, and (b) representation of a confidence region for the MLE: that is what distinguishes a statistical analysis from a purely numerical fitting problem.
– whuber
Commented Dec 16, 2023 at 16:49
• Thanks for the comments, I'll update this answer later today Commented Dec 16, 2023 at 17:15

I think you’ll need to assign a distribution to the observed distance (or the square of the distance) that has a mean and variance dependent on the true distance. I don’t know the characteristics of the devices you’re using and offer the following approach only as a possible structure to follow.

I’m assuming from your description that from each receiver you only obtain a distance and NOT a direction.

Suppose the distribution of the square of the observed distance follows a Gamma distribution with the mean equal to the square of the true distance and the variance is proportional to the square of the true distance. (I’m assuming that smaller distances are more precisely estimated than large differences.)

If the true squared distance is $$\mu^2$$ and the variance is $$\mu^2 \sigma^2$$, then the associated Gamma distribution parameters are $$\alpha=\mu^2/\sigma^2$$ and $$\beta=\sigma^2$$.

Below is Mathematica code for generating data and finding the maximum likelihood estimates.

(* Transmitter location *)
a = {0.75, 0.8};

r = {{0.25, 0.25}, {0.25, 0.75}, {0.75, 0.25}, {0.75, 0.75}};
n = Length[r] (* Number of receivers *)
trueDist2 = Table[Norm[r[[i]] - a]^2, {i, n}] (* True squared distances to transmitter *)
trueσ = 0.05;  (* Error parameter *)

(* Sample of observations of squared distance *)
SeedRandom[12345];
d2 = Table[RandomVariate[GammaDistribution[trueDist2[[i]]/trueσ^2, trueσ^2]], {i, n}]
(* {0.5525,0.2525,0.3025,0.0025} *)

(* Log of the likelihood *)
logL = Sum[LogLikelihood[GammaDistribution[((r[[i, 1]] - μx)^2 + (r[[i, 2]] - μy)^2)/σ^2, σ^2], {d2[[i]]}],
{i, n}];

(* Maximum likelihood estimates *)
mle = FindMaximum[{logL, σ > 0}, {{μx, a[[1]]}, {μy, a[[2]]}, {σ, trueσ}}, MaxIterations -> 5000]
(* {11.571356964839643,{μx -> 0.7146660584863823, μy -> 0.7750478956280523, σ -> 0.04110551720224718}} *)

(* Estimate of covariance matrix *)
(cov = -Inverse[D[logL, {{μx, μy, σ}, 2}] /. mle[[2]]]) // MatrixForm
se = Diagonal[cov]^0.5  (* Standard errors of estimates of μx, μy, and σ *)
(* {0.01457950738463684, 0.01626917105787927, 0.01464057953398463} *)


If it helps, here is equivalent R code:

# Transmitter location
a <- c(0.75, 0.8)

n <- 4  # Number of recieviers
r <- matrix(c(0.25, 0.25, 0.25, 0.75, 0.75, 0.25, 0.75, 0.75), nrow=4, byrow=TRUE)

# True squared distances to transmitter
trueDist2 = (r[,1] - a[1])^2 + (r[,2] - a[2])^2
# [1] 0.5525 0.2525 0.3025 0.0025

# Error parameter associated with measurement of squared distance
truesigma <- 0.05

# Sample of observations of squared distance
set.seed(12345)
d2 <- rgamma(n, shape=trueDist2/truesigma^2, scale=truesigma^2)
d2 <- c(0.454823, 0.24621, 0.286231, 0.00226926)  # Squared distances to match Mathematica example

# Define log likelihood function
logL <- function(parms, r, d2) {
mx <- parms[1]
my <- parms[2]
sigma <- parms[3]
temp <- ((r[, 1] - mx)^2 + (r[, 2] - my)^2)/sigma^2
-sum(d2)/sigma^2 - sum(lgamma(temp)) - sum(log(sigma^2)*temp) + sum(log(d2)*(-1 + temp))
}

# Maximize log likelihood
result <- optim(c(a[1], a[2], truesigma), logL, r=r, d2=d2,
control=list(fnscale=-1), hessian=TRUE)
result

# Covariance matrix and standard errors for parameters
cov <- -solve(resulthessian) se <- diag(cov)^0.5 # {0.01457950738463684, 0.01626917105787927, 0.01464057953398463} *)  Determining the equations for the MLE can be done with standard calculus methods. To facilitate our analysis we can define the functions: $$r_i \equiv r_i (x_b, y_b) \equiv \sqrt{(x_i - x_b)^2+(y_i - y_b)^2}.$$ The log-likelihood function for your model is: $$\ell(x_b, y_b, \sigma) = \text{const} - n \log(\sigma) - \frac{1}{2 \sigma^2} \sum_i (r_i - d_i)^2.$$ The partial derivatives of this log-likelihood are: \begin{align} \frac{\partial \text{SE}}{\partial x_b} (x_b, y_b, \sigma) &= \sum_i \frac{2(x_i - x_b) (r_i - d_i)^2}{r_i}, \\[12pt] \frac{\partial \text{SE}}{\partial y_b} (x_b, y_b, \sigma) &= \sum_i \frac{2(y_i - y_b) (r_i - d_i)^2}{r_i}, \\[12pt] \frac{\partial \text{SE}}{\partial \sigma} (x_b, y_b, \sigma) &= - \frac{n}{\sigma} + \frac{1}{\sigma^3} \sum_i (r_i - d_i)^2. \\[12pt] \end{align} Setting the partial derivatives equal to zero gives the score equation for the MLE: \begin{align} 0 &= \sum_i \frac{2(x_i - x_b) (\hat{r}_i - d_i)^2}{\hat{r}_i}, \\[12pt] 0 &= \sum_i \frac{(y_i - y_b) (\hat{r}_i - d_i)^2}{\hat{r}_i}, \\[12pt] \hat{\sigma}^2 &= \frac{1}{n} \sum_i (\hat{r}_i - d_i)^2. \\[12pt] \end{align} There is no closed-form solution to these equations so you will need to solve them numerically (e.g., with Newton-Raphson iteration). (Note that the MLE for the parameter $$\sigma$$ is typically biased and is usually replaced with an alternative bias-mitigated estimator using an adjustment like Bessel's correction.) • Should it be made explicit that one is assuming that the error is at least approximately normal with a constant variance? Even though the OP doesn't mention it, it would seem likely or at least possible that larger distances are associated with a larger error. – JimB Commented Dec 16, 2023 at 17:06 • @JimB: OP already states this explicitly in the question (\epsilon_i \sim \mathcal{N}(0,\,\sigma^{2})\$).
– Ben
Commented Dec 16, 2023 at 20:00
• Sorry. My feeble excuse is that I didn't re-read the question after @User1865345 added that in.
– JimB
Commented Dec 16, 2023 at 20:58
• @JimB: All good. No apology required.
– Ben
Commented Dec 16, 2023 at 22:02