What is the Maximum Likelihood Estimator for my nonlinear problem?

Question

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?

Answer 1

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.

Answer 2

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"))
    # add legend
    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)

Answer 3

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};

(* Receiver locations *)
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)

# Receiver locations
  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(result$hessian)
  se <- diag(cov)^0.5 
  # {0.01457950738463684, 0.01626917105787927, 0.01464057953398463} *)     

Answer 4

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.)