I'm doing some research into the Cramer-Rao bound for time of arrival localization and have come across a rather strange result: the FIM is singular, but there exists an unbiased estimator. My supervisor insists I'm wrong (I'm sure I must be), but I can't seem to figure out what I'm doing wrong. Can someone please help guide me in the right direction? Sorry for the long development of the problem below.
A transmitter $T$ is at an unknown location $\theta=(x, y)^T$. There are $s$ sensors in the region $\Omega$. The $j^\text{th}$ sensor positioned at $\theta_j=(x_j, y_j)^T$ measures a range to $T$ that is a random variable dependent on the actual distance of $T$ from the sensor ($d_j$), as shown below:
The range measurement vector for all anchors is therefore: $$ r = [d_1 + \varepsilon_1, \ldots, d_j + \varepsilon_j]^T = d(\theta)+\varepsilon=\|\theta-\theta_j\|_2 + \varepsilon $$
where the error is assumed to be normally distributed $\epsilon \sim N\left(0, \Sigma_{s_{\mathrm{x}}s} \right)$.
For a multivariate Gaussian distribution, the $(m,n)$ element of the FIM is:
$$ FIM_{m,n}=\frac{\partial d^T}{\partial \theta_m} \Sigma^{-1} \frac{\partial d}{\partial \theta_n}+\frac{1}{2} \operatorname{TR}\left(\Sigma^{-1} \frac{\partial \Sigma}{\partial \theta_m} \Sigma^{-1} \frac{\partial \Sigma}{\partial \theta_n}\right) $$
Assuming independent noise at the sensor measurement and a constant variance for the noise at each sensor, the covariance matrix is diagonal, i.e. $\Sigma=\operatorname{diag}(\sigma^2_1, \sigma^2_2,\ldots, \sigma^2_s)$. Thus, the FIM can be expanded for $\theta=(x, y)^T$ as:
\begin{equation} \left[\begin{array}{cc} FIM_{1,1} & FIM_{1,2} \\ FIM_{2,1} & FIM_{2,2} \end{array}\right] \end{equation}
\begin{align} \label{FIM11} FIM_{1,1} &= \sum_{j=1}^s \frac{1}{2} \left(\sigma_j^2\right)^{-2} \left(\frac{\partial \sigma_j^2}{\partial x} \right)^2 + \left(\sigma_j^2\right)^{-1}\left(\frac{\partial d_j}{\partial x}\right)^2 \\ & = \sum_{j=1}^s \left(\sigma_j^2\right)^{-1} \cos^2 (\phi_j) \end{align}
\begin{align} \label{FIM12} FIM_{1,2}=FIM_{2,1}&=\sum_{j=1}^{s} \frac{1}{2} \left(\sigma_j^2\right)^{-2} \frac{\partial \sigma_j^2}{\partial x} \frac{\partial \sigma_j^2}{\partial y} + \left(\sigma_j^2 \right)^{-1} \frac{\partial d_j}{\partial x} \frac{\partial d_j}{\partial y}\\ &=\sum_{j=1}^s \left(\sigma_j^2 \right)^{-1} \cos(\phi_j) \sin(\phi_j) \end{align}
\begin{align} \label{FIM22} FIM_{2,2} &= \sum_{j=1}^s \frac{1}{2}\left(\sigma_j^2 \right)^{-2} \left(\frac{\partial \sigma_j^2}{\partial y}\right)^2 + \left(\sigma_j^2 \right)^{-1} \left(\frac{\partial d_j}{\partial y}\right)^2 \\ &=\sum_{j=1}^s \left(\sigma_j^2 \right)^{-1} \sin(\phi_j)^2 \end{align}
where $\phi_j$ is the angle between the $j^\text{th}$ sensor and T, since
\begin{align} \frac{\partial d_{j}}{\partial x}&=\frac{\partial}{\partial x} \left[(x-x_j)^{2}+(y-y_j)^{2}\right]^{\frac{1}{2}}=(x-x_j)d_{j}^{-1}\\ &=cos(\phi_j) \end{align}
And similarly $\frac{\partial d_j}{\partial y}=sin(\phi_j)$.
Consider three sensors situated at $(-1,0),$ $(0,0)$ and $(1,0),$ and $T$ at $(2,0).$ The FIM is clearly singular as $\phi_j=0 \, \forall j$ and so the CRB does not exist. The negative log-likelihood function (NLL) can be written as:
\begin{align} -\mathcal{L}&=-\log\left(\prod_{j=1}^{s} P_{r}\left(r_{j} ; d_{j}\right)\right) = -\log\left(\prod_{j=1}^{s} \frac{1}{\sqrt{2 \pi \sigma_{j}^{2}}} e^{-\frac{\left(r_{j}-d_{j}\right)^{2}}{2 \sigma_{j}^{2}}}\right)\\ &= \sum_{j=1}^{s}\left[\frac{1}{2}\log (2 \pi)+\frac{1}{2} \log \left(\sigma^{2}_{j}\right)+\frac{(r-d_{j})^{2}}{2 \sigma^{2}_{j}}\right] \end{align}
I've done an MLE simulation in MATLAB for this.
The circles for two different range measurements vectors are shown below, overlayed with the negative log-likelihood (NLL) function for the measurement.
There are two situations that happen with the range circles:
When 2 or more circles intersect, they do so at two locations (the two bluish regions in the left plot) which have reflection symmetry along the $y=0$ line, and so the NLL function has a minima at these two points.
When no circles intersect, the NLL function has a minima on the $y=0$ line (think of this in a least-squares sense: when the circles don't intersect, the least squares estimate is where they are closest, which is the bluish region in the right plot).
Therefore, the $y$-coordinate MLE (i.e. estimator that minimizes NLL function) is unbiased and has finite variance.
This is the output of 1500 MLE estimates I computed in a simulation:
I can't seem to figure out where I've gone wrong.
EDIT: I asked this question in math.stackexchange, but I figure this is probably a better place to ask it since it's PhD research that I'm doing. Also as background, I'm an engineer but part of my supervisory team looking at this part of my work is composed of mathematicians.