# An Approximate MSE Expression for Maximum Likelihood and Other Implicitly
Defined Estimators of Non-Random Parameters
(extended version)

Erdal Mehmetcik<EMAIL_ADDRESS>Umut Orguner<EMAIL_ADDRESS>Çağatay Candan<EMAIL_ADDRESS>Department of Electrical and Electronics
Engineering, Middle East Technical University (METU), 06800, Ankara, Turkey.

###### Abstract

An approximate mean square error (MSE) expression for the performance analysis
of implicitly defined estimators of non-random parameters is proposed. An
implicitly defined estimator (IDE) declares the minimizer/maximizer of a
selected cost/reward function as the parameter estimate. The maximum
likelihood (ML) and the least squares estimators are among the well known
examples of this class. In this paper, an exact MSE expression for implicitly
defined estimators with a symmetric and unimodal objective function is given.
It is shown that the expression reduces to the Cramer-Rao lower bound (CRLB)
and misspecified CRLB in the large sample size regime for ML and misspecified
ML estimation, respectively. The expression is shown to yield the Ziv-Zakai
bound (without the valley filling function) for the maximum a posteriori (MAP)
estimator when it is used in a Bayesian setting, that is, when an a-priori
distribution is assigned to the unknown parameter. In addition, extension of
the suggested expression to the case of nuisance parameters is studied and
some approximations are given to ease the computations for this case.
Numerical results indicate that the suggested MSE expression not only predicts
the estimator performance in the asymptotic region; but it is also applicable
for the threshold region analysis, even for IDEs whose objective functions do
not satisfy the symmetry and unimodality assumptions. Advantages of the
suggested MSE expression are its conceptual simplicity and its relatively
straightforward numerical calculation due to the reduction of the estimation
problem to a binary hypothesis testing problem, similar to the usage of Ziv-
Zakai bounds in random parameter estimation problems.

###### keywords:

Parameter Estimation , Implicitly Defined Estimator , Maximum Likelihood ,
Misspecified Maximum Likelihood , Performance Prediction , Cramer-Rao Lower
Bound , Ziv-Zakai Bound.

††journal: Signal Processing

##  1 Introduction

The topic of parameter estimation can be divided into two classes, namely the
estimation of random and non-random (deterministic) parameters. The random
parameter estimation (Bayesian estimation) assumes that the parameter of
interest is a random variable with an a-priori distribution and the
observations on a realization of the unknown parameter are obtained according
to a known probabilistic mapping. Under this setting, the optimal estimator
that minimizes the risk, say mean square error (MSE) or mean absolute error,
is a functional of the posterior density of the parameter. For instance, the
optimal estimator minimizing the MSE is the mean value of the parameter with
respect to the posterior density [1]. In general, the posterior density
calculation is the key step for the Bayesian formulation. Unfortunately, a
closed form expression for the posterior density (and its moments) which does
not involve integration, differentiation and limit operations is rarely
available. In many problems, one has to resort to the Monte Carlo methods or
approximate inference techniques for an inexact realization of the optimal
Bayesian estimator. In such problems, the estimator success is typically
evaluated by comparisons with the performance bounds. Bayesian performance
bounds have a vast literature [2]. Typically, these bounds do not impose any
constraints on the estimator. For instance the Bayesian Cramer Rao lower bound
(CRLB) [2, 3], Weiss Weinstein bound (WWB) [4], Bayesian Bhattacharya bound
[3] are derived using the covariance inequality principle (hence, sometimes
referred as covariance bounds) and applicable in general to any type of
estimators. Another main class is the Ziv-Zakai bound (ZZB) [5, 6, 7] type
bounds which are derived by converting the estimation problem into a binary
detection problem. Bayesian CRLB is one of most fundamental bounds and
provides the achievable MSE in the asymptotic region which is the high signal-
to-noise ratio (SNR) region. However, it suffers from the threshold effect [5,
6, 7], meaning that it provides unachievable (optimistic) lower bounds at
medium or low SNR values. ZZB and WWB are among the tightest Bayesian bounds
in all regions of operation [2, 7]. Bayesian bounds continue to be an active
research area. Recently, Bayesian bounds for estimating periodic parameters
(e.g., phase) have been developed [8, 9, 10].

Non-random parameter estimation involves some challenges unique to this
setting. In this setting, an estimator can be improved for a specific value of
the parameter at the expense of performance for other parameter values [1,
11]. Since there is no a-priori distribution associated to the parameter of
interest, it is not possible to balance the performance gains and losses for
different parameter values as in the Bayesian setting. For example, the
estimator ignoring the measurements and producing a constant value, say
$\alpha$, as the estimate has no error if the unknown parameter is indeed
$\alpha$; but, suffers from performance losses at all other parameter values.
The development of lower bounds for the non-random parameter estimation also
suffers from similar inherent admissibility problems. To overcome these
problems, the estimators in this setting are typically restricted to the class
of unbiased estimators and examined under the title of minimum variance
unbiased estimators [1].

The performance bounds for the non-random parameter estimation are also
developed for a specific class of estimators. For example CRLB (for non-random
parameters), Hammersley-Chapman-Robbins Bound (HCRB) [12], Barankin Bound [13,
14] require the estimator to be unbiased in an open neighborhood of a point,
over a set of two-points and over a set of many-points, respectively (also see
[15]). In [16], a general bound form for unbiased estimators is given and it
is shown that CRLB, HCRB and BB can be derived by a proper choice of the
kernel function of their integral transform. Note that the unbiasedness
condition may not be practical or may be difficult to satisfy, especially for
the parameters with a finite support; since the estimation error approaches a
one-sided distribution at the edges of the parameter space in such cases [5].
Although it has been shown that for some problems with periodic parameters
[17, 18, 19], the problem with the one-sided error distribution at the edges
may vanish; uniformly unbiased estimators do not exist in these cases either
[20]. Furthermore, the unbiasedness condition may not even be desirable in
some problems. It is known that there exist realizable biased-estimators for
some problems whose MSE is lower than the Cramer-Rao bound for unbiased
estimators [21]. Perhaps, the most important aspect of unbiasedness condition
is in relation with the maximum likelihood estimator. It is well known that
the maximum likelihood estimator is unbiased and efficient in the large sample
size regime, under fairly general conditions, providing a basis for the
theoretical and practical adoption of the unbiasedness condition [22, 23].

The main problem considered in this paper is the performance prediction of
implicitly defined estimators (IDEs) of non-random parameters. IDEs are
estimators which produce an estimate by maximizing an objective function of
the measurements over the parameter set under consideration. The maximum
likelihood (ML) estimator, least squares estimators are some well known
examples.

In this study, we present an approximate MSE expression for IDEs of non-random
parameters that

  * 1.

gives the true MSE when the objective function of the IDE is symmetric and
unimodal,

  * 2.

reduces to the CRLB in the large sample size regime for ML estimation,

  * 3.

reduces to the misspecified CRLB (MCRLB) [24, 25] in the large sample size
regime for misspecified ML estimation [25, 26],

  * 4.

reduces to the ZZB when an a-priori distribution is assigned to the parameter
of interest for maximum a posteriori (MAP) estimation.

There are already some approximate MSE expressions available in the literature
for IDEs of non-random parameters. For instance, [27] provides formulas for
MSE and bias of IDEs using Taylor series expansion of the cost function along
with some approximate expressions for certain expectations and derivatives.
The study in [28] also uses Taylor expansion approach and derives different
approximations in the scalar parameter case. Both of these approaches are
based on the Taylor series expansion around the true parameter value and
provide simple MSE expressions; but, do not take into account the gross errors
which becomes the significant factor as the estimator nonlinearity increases
and/or SNR is decreased below the threshold SNR.

The majority of work on the estimator performance prediction focus on the
performance of the ML estimator. The ML estimator is known to be
asymptotically efficient (performance approaching CRLB at large sample size)
under some regularity conditions [3, 22, 23]. The method of interval errors
(MIE) is a celebrated method that was proposed by Van Trees [3] to assess the
performance of the ML estimator in the threshold region. This method depends
on a careful selection of intervals in the parameter space and the calculation
of their probabilities. Different approximations have been proposed to
approximate the probabilities [29, 30, 31]. The MSE expression proposed in the
present study can be interpreted as a more principled version of the method of
interval errors where the need for the interval selection and the gross error
probability calculation or approximation is not required.

Notation: Throughout the paper lower and uppercase letters denote scalars,
e.g., $a$, $A$. Bold lowercase letters denote vectors, e.g., $\mathbf{a}$.
Bold uppercase letters denote matrices, e.g., $\mathbf{A}$. The $i$th element
of the vector $\mathbf{a}$ is denoted by $[\mathbf{a}]_{i}$. The $i,j$th
element of the matrix $\mathbf{A}$ is denoted by $[\mathbf{A}]_{i,j}$.
$\Re\\{\cdot\\}$ denotes the real part of the complex argument.

##  2 Problem Definition

We consider the estimation of the non-random real-valued vector
$\boldsymbol{\theta}\triangleq[\theta_{1}\,\,\theta_{2}\,\cdots\,\,\theta_{J}]^{\mathrm{T}}{}$
from the measurements
$\mathbf{x}\triangleq[x_{0}\,\,x_{1}\,\cdots\,\,x_{N-1}]^{\mathrm{T}}{}\in\mathbb{C}^{N}$
distributed according to $f(\mathbf{x};\boldsymbol{\bar{\theta}})$ where
$\boldsymbol{\bar{\theta}}\triangleq[\bar{\theta}_{1}\,\,\bar{\theta}_{2}\,\cdots\,\,\bar{\theta}_{J}]^{\mathrm{T}}{}$
denotes the true value of $\boldsymbol{\theta}$. An implicitly defined
estimator (IDE) generates an estimate
$\boldsymbol{\hat{\theta}}\triangleq[\hat{\theta}_{1}\,\,\hat{\theta}_{2}\,\cdots\,\,\hat{\theta}_{J}]^{\mathrm{T}}{}$
by maximizing an objective function $\mathcal{L}(\cdot,\cdot)$ of the
measurements and the parameters as shown below:

$\displaystyle\boldsymbol{\hat{\theta}}\triangleq\arg\max_{\boldsymbol{\theta}}\mathcal{L}(\mathbf{x};\boldsymbol{\theta}).$
(1)

The most well-known example of IDEs is the ML estimator where the objective
function $\mathcal{L}(\cdot,\cdot)$ is the likelihood function
$f(\mathbf{x};\boldsymbol{\theta})$. Other examples of IDEs are M-estimators
and (nonlinear) least square estimators. We see that the estimate
$\boldsymbol{\hat{\theta}}$ given by (1) is determined by the measurements
implicitly, hence the name _implicitly defined estimator_.

In this study we are interested in the performance of IDEs and we give an
expression for the (diagonal elements of the) MSE matrix of the estimate
$\boldsymbol{\hat{\theta}}$ which is defined as

$\displaystyle\operatorname{MSE}(\boldsymbol{\bar{\theta}})\triangleq$
$\displaystyle\mathbb{E}\big{[}(\boldsymbol{\hat{\theta}}-\boldsymbol{\bar{\theta}})(\boldsymbol{\hat{\theta}}-\boldsymbol{\bar{\theta}})^{\mathrm{T}}\big{]},$
(2)

Here it should be mentioned that the methodology presented in the current work
can be straightforwardly extended to the other moments of the estimation error
$\boldsymbol{\hat{\theta}}-\boldsymbol{\bar{\theta}}$.

Except for few cases like ML estimation for Gaussian likelihoods with linear
models, the estimate $\boldsymbol{\hat{\theta}}$ in (1) cannot be analytically
expressed in terms of the measurements $\mathbf{x}$, i.e., one cannot find a
closed form expression for the function $\mathbf{h}(\cdot)$ such that
$\boldsymbol{\hat{\theta}}=\mathbf{h}(\mathbf{x})$. As a consequence the
determination, evaluation and comparison of performance (say, in terms of MSE)
of an IDE usually involves extensive Monte Carlo studies and/or problem
specific approximations. In this work we first give an MSE expression which is
exact for an IDE of a scalar parameter whose objective function is both
symmetric (around the estimate) and unimodal in Section 3. Since the symmetry
and unimodality conditions are typically satisfied by the objective functions
of IDEs in the asymptotic or small error region, as further examined in
Section 4; we suggest to use the MSE expression to study the performance of
IDEs in the small error and threshold regions. We refrain from calling the
suggested MSE expression as a bound due to the lack of performance guarantees
in the non-asymptotic regions. The suggested expression can be considered to
be in the same league with the MIE [3] which lacks a performance guarantee in
all regions including the asymptotic region. Such expressions are also called
approximate bounds in some studies [31, 32]. Yet, our main goal in this study
is to develop an MSE expression similar to ZZB, which is known to be a tight
random parameter estimation bound in the threshold and asymptotic regions, for
non-random parameters.

##  3 Case of a Scalar Parameter with Symmetric and Unimodal Objective
Functions

In this section we are going to restrict ourselves to a scalar unknown
parameter $\theta\in\mathbb{R}$ (i.e., $J=1$) and provide a predicted MSE
expression which is equal to the true MSE for an IDE whose objective function
satisfies symmetry and unimodality assumptions. Our main results are given in
the following theorem and its corollary.

Figure 1: Illustration of the fact that the events
$(\hat{\theta}-\theta)\geq\epsilon$ and $(\hat{\theta}-\theta)\leq-\epsilon$
are equivalent to the events
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$
and
$\mathcal{L}(\mathbf{x};\theta-2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$,
respectively, when the assumptions of Theorem 1 hold.

######  Theorem 1.

Consider the IDE given as

$\displaystyle\hat{\theta}\triangleq\arg\max_{\theta}\mathcal{L}(\mathbf{x};\theta).$
(3)

Let the objective function $\mathcal{L}(\mathbf{x};\cdot)$ satisfy the
following conditions for all $\mathbf{x}\in\mathbb{C}^{N}$.

  1. 1.

$\mathcal{L}(\mathbf{x};\hat{\theta}+h)=\mathcal{L}(\mathbf{x};\hat{\theta}-h)$
for all $h\in\mathbb{R}$, i.e., the objective function is symmetric around its
peak.

  2. 2.

$\mathcal{L}(\mathbf{x};\theta)$ is strictly-increasing (strictly-decreasing)
for $\theta<\hat{\theta}$ ($\theta>\hat{\theta}$).

Define the true estimator statistic $V_{\hat{\theta}}(\theta)$ as

$\displaystyle V_{\hat{\theta}}(\theta)\triangleq
E[(\hat{\theta}-\theta)^{2}],$ (4)

where $\theta$ is an arbitrary fixed parameter value. Then,

$\displaystyle V_{\hat{\theta}}(\theta)$
$\displaystyle=\widehat{V}_{\hat{\theta}}(\theta),$ (5)

where the predicted statistic $\widehat{V}_{\hat{\theta}}(\theta)$ is defined
as

$\displaystyle\widehat{V}_{\hat{\theta}}(\theta)$ $\displaystyle\triangleq
2\int_{-\infty}^{\infty}\\!|\epsilon|P\left(\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)\\!\\!\;\mathrm{d}{}\epsilon.$
(6)

Proof: A proof is presented in A. The main idea of the proof is to show that
the events $(\hat{\theta}-\theta)\geq\epsilon$ and
$(\hat{\theta}-\theta)\leq-\epsilon$ are equivalent to the events
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$
and
$\mathcal{L}(\mathbf{x};\theta-2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$,
respectively, when the assumptions of the theorem hold. An illustration of
this equivalence is given in Figure 1. $\blacksquare$

The following remark applies Theorem 1 to find the true MSE of the estimator
$\hat{\theta}$.

######  Remark 1 (MSE of IDE).

The true MSE of the IDE $\hat{\theta}$ in Theorem 1 is given as

$\displaystyle\operatorname{MSE}(\bar{\theta})$
$\displaystyle=\widehat{\operatorname{MSE}}(\bar{\theta})$ (7)

where the predicted MSE, denoted as
$\widehat{\operatorname{MSE}}(\bar{\theta})$, is defined as

$\displaystyle\widehat{\operatorname{MSE}}(\bar{\theta})$
$\displaystyle\triangleq
2\int_{-\infty}^{\infty}\\!|\epsilon|P\left(\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})\right)\\!\\!\;\mathrm{d}{}\epsilon.$
(8)

Proof: The proof is trivial by realizing that
$\operatorname{MSE}(\bar{\theta})=V_{\hat{\theta}}(\bar{\theta})$ and
$\widehat{\operatorname{MSE}}(\bar{\theta})=\widehat{V}_{\hat{\theta}}(\bar{\theta})$.
$\blacksquare$

In the special case of a parameter $\theta$ with finite support, e.g.,
$\theta\in[\theta_{\min},\,\theta_{\max}]$, the estimation error
$\hat{\theta}-\bar{\theta}$ is restricted to the interval
$[\theta_{\min}-\bar{\theta},\,\theta_{\max}-\bar{\theta}]$ and the predicted
MSE becomes

$\displaystyle\widehat{\operatorname{MSE}}(\bar{\theta})=\,2\int_{\frac{\theta_{\min}-\bar{\theta}}{2}}^{\frac{\theta_{\max}-\bar{\theta}}{2}}|\epsilon|P\left(\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})\right)\;\mathrm{d}{}\epsilon.$
(9)

We can interpret the MSE expression (8) intuitively as follows. When the
probability
$P\left(\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})\right)$
is large for large values of $|\epsilon|$, then it is probable for the IDE
$\hat{\theta}$ to make gross errors, resulting in a large MSE. On the other
hand, if this probability is small for large values of $|\epsilon|$, the
contribution of gross errors in the MSE becomes negligible, resulting in a
small MSE. Consequently, IDEs with a small MSE would have the probability
$P\left(\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})\right)$
(thought of as a function of $\epsilon$) highly concentrated in a small
neighborhood of $\epsilon=0$ and _quickly_ vanishing elsewhere. More
specifically, a sufficient and necessary condition for existence of the
integral in the MSE expression (8) is
$P\left(\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})\right)=o(1/|\epsilon|^{2})$,
i.e., the probability
$P\left(\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})\right)$
decaying strictly faster than $1/|\epsilon|^{2}$ as
$|\epsilon|\rightarrow\infty$. The integral in the MSE expression (9), on the
other hand, always exists.

We can put the expression $\widehat{\operatorname{MSE}}(\bar{\theta})$ to a
test by considering the optimal but infeasible estimator
$\hat{\theta}=\bar{\theta}$. This estimator can be formulated as an IDE using
the objective function
$\mathcal{L}(\mathbf{x};\theta)\triangleq-(\theta-\bar{\theta})^{2}$. Since
the objective function $\mathcal{L}(\mathbf{x};\cdot)$ does not depend on the
measurements $\mathbf{x}$, we see that the probability of the deterministic
event
$\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})$
is given as

$\displaystyle
P\left(\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})\right)=\begin{cases}1,&\epsilon=0\\\
0,&\text{otherwise}\end{cases}.$ (10)

When (10) is substituted into (8), we have
$\widehat{\operatorname{MSE}}(\bar{\theta})=0$, which is the true MSE.
Similarly, for the feasible (but biased) version of this estimator
$\hat{\theta}=\theta_{0}$, where $\theta_{0}\in\mathbb{R}$, the objective
function is $\mathcal{L}(\mathbf{x};\theta)\triangleq-(\theta-\theta_{0})^{2}$
and the corresponding probability becomes

$P\left(\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})\right)=\begin{cases}1,&-(\bar{\theta}+2\epsilon-\theta_{0})^{2}\geq-(\bar{\theta}-\theta_{0})^{2}\\\
0,&\text{otherwise}\end{cases}=\begin{cases}1,&0\geq\epsilon\geq-\bar{\theta}+\theta_{0}\\\
1,&-\bar{\theta}+\theta_{0}\geq\epsilon\geq 0\\\
0,&\text{otherwise}\end{cases},$ (11)

which yields
$\widehat{\operatorname{MSE}}(\bar{\theta})=(\bar{\theta}-\theta_{0})^{2}$
when substituted into (8) for both $\theta_{0}\geq\bar{\theta}$ and
$\theta_{0}\leq\bar{\theta}$. This also is the true MSE. The following
corollary applies the result in Theorem 1 to ML estimation.

######  Corollary 1 (MSE of ML Estimator).

If the likelihood function $f(\mathbf{x};\cdot)$ satisfies the conditions in
Theorem 1, then the true MSE of the ML estimate $\hat{\theta}$ is given as

$\displaystyle\operatorname{MSE}_{\mathrm{ML}}(\bar{\theta})=\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})$
(12)

where the predicted MSE, denoted as
$\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})$, is defined as

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})\triangleq$
$\displaystyle\,2\int_{-\infty}^{\infty}|\epsilon|P\bigg{(}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
1\bigg{)}\;\mathrm{d}{}\epsilon.$ (13)

Proof: Since we have $\mathcal{L}(\mathbf{x};\theta)\triangleq
f(\mathbf{x};\theta)$, we can write

$\displaystyle
P\left(\mathcal{L}(\mathbf{x};\bar{\theta}+2\epsilon)\geq\mathcal{L}(\mathbf{x};\bar{\theta})\right)$
$\displaystyle=P\left(f(\mathbf{x};\bar{\theta}+2\epsilon)\geq
f(\mathbf{x};\bar{\theta})\right)=\int
I\left(f(\mathbf{x};\bar{\theta}+2\epsilon)\geq
f(\mathbf{x};\bar{\theta})\right)f(\mathbf{x};\bar{\theta})\;\mathrm{d}{}\mathbf{x}$
(14a) $\displaystyle=\int_{f(\mathbf{x};\bar{\theta})\neq
0}I\left(f(\mathbf{x};\bar{\theta}+2\epsilon)\geq
f(\mathbf{x};\bar{\theta})\right)f(\mathbf{x};\bar{\theta})\;\mathrm{d}{}\mathbf{x}=\int_{f(\mathbf{x};\bar{\theta})\neq
0}I\bigg{(}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
1\bigg{)}f(\mathbf{x};\bar{\theta})\;\mathrm{d}{}\mathbf{x}$ (14b)
$\displaystyle=P\bigg{(}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
1\bigg{)},$ (14c)

where $I(\cdot)$ denotes the indicator function for event arguments.
Substituting the last probability into the integral in (8) completes the
proof. $\blacksquare$

Note that the expression (13) connects the MSE of the ML estimate to the error
probability of a likelihood ratio test. This connection between estimation and
detection theory is further explored in Section 4.3 in relation with the ZZB.

Theorem 1 and its corollary provide compact expressions to evaluate the MSE of
an implicitly defined estimator exactly even when there is no explicit
analytical expression connecting the estimate $\hat{\theta}$ to the
measurements $\mathbf{x}$. However, it has some limitations imposed by the
assumptions required for its validity. In fact, almost all practical
estimation problems violate one of the assumptions of symmetry, unimodality
and infinite support (of the parameter $\theta$). For these problems it is
certainly possible to have
$\operatorname{MSE}(\bar{\theta})\neq\widehat{\operatorname{MSE}}(\bar{\theta})$.
Hence, in general the proposed expressions in (8) and (13) can only serve as
approximate MSE performance prediction tools. Keeping this fact in mind, we
show several relations between the suggested MSE expression and well-known
bounds in Section 4.

##  4 Relationship to Performance Bounds

In this section we present the relationship of the suggested MSE expression
(13) to some well-known performance bounds.

###  4.1 Relationship to CRLB

In this section, we consider the ML estimation for a scalar parameter
$\theta\in\mathbb{R}$. In order to use the large sample asymptotic results for
the ML estimate $\hat{\theta}$, we are going to assume that the elements
$x_{n}$, $n=0,\ldots,N-1$, of the measurement vector $\mathbf{x}$ are
independent and identically distributed as $x_{n}\sim f(x_{n};\bar{\theta})$
111In order not to incorporate additional notation, we will keep here the
individual measurements $x_{n}\in\mathbb{C}$ as scalars but the same results
can be obtained for the case when $x_{n}$ is a vector.. The likelihood for the
measurement vector $\mathbf{x}$ is then given as

$\displaystyle f(\mathbf{x};\theta)=\prod_{n=0}^{N-1}f(x_{n};\theta).$ (15)

The relationship of the suggested MSE expression
$\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})$ to CRLB is given in
the following proposition.

######  Proposition 1.

  * Assume that

  * A0

The parameter $\theta$ has finite support, i.e.,
$\theta\in[\theta_{\mathrm{min}},\theta_{\mathrm{max}}]$, and the true
parameter value $\bar{\theta}$ satisfies
$\bar{\theta}\in(\theta_{\mathrm{min}},\theta_{\mathrm{max}})$.

  * A1

The first three derivatives of $\ln f(x;\theta)$ with respect to $\theta$
exist for all $\theta$ and are continuous with respect to $\theta$.

  * A2

For every $\theta$, the functions $|\frac{\partial^{i}}{\partial\theta^{i}}\ln
f(x;\theta)|$, $i=0,1,2,3$, are dominated by functions $b_{i}(x)$,
$i=0,1,2,3$, which all have finite variance.

  * A3

The KL divergence $D(f(x;\bar{\theta})||f(x;\theta))$, where

$\displaystyle D(f(x)||g(x))\triangleq\int
f(x)\ln\frac{f(x)}{g(x)}\;\mathrm{d}{}x,$ (16)

has a unique minimum with respect to $\theta$ at $\theta=\bar{\theta}$.

  * A4

The expectation $\mathbb{E}\big{[}\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(x;\bar{\theta})\big{]}$ is non-zero.

The assumption A0 is sufficient (but not necessary) for the convergence of the
integral in (13). Under the assumptions A1-A3, it can be shown that (See [22,
Theorem 2.1]) the ML estimate $\hat{\theta}\triangleq\arg\max_{\theta}\ln
f(\mathbf{x};\theta)$ is consistent, i.e.,
$\hat{\theta}\overset{\text{a.s.}}{\rightarrow}\bar{\theta}$ as
$N\rightarrow\infty$ . Then, we have

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})\rightarrow\operatorname{MSE}_{\mathrm{ML}}(\bar{\theta})$
(17)

as $N\rightarrow\infty$, i.e., the finite support version of the MSE
expression $\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})$ in (13)
converges to the true large sample asymptotic MSE of the ML estimate
$\hat{\theta}$ as $N\rightarrow\infty$. If the ML estimate is also
asymptotically efficient, then we have

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})\rightarrow\mathcal{C}(\bar{\theta}),$
(18)

as $N\rightarrow\infty$ where
$\mathcal{C}(\bar{\theta})\triangleq\mathcal{I}^{-1}(\bar{\theta})$ with
$\mathcal{C}(\bar{\theta})$ and $\mathcal{I}(\bar{\theta})$ denoting the CRLB
and the Fisher information matrix, respectively, at the true parameter value
$\bar{\theta}$.

Proof: The proof is given in B. $\blacksquare$

###  4.2 Relationship to MCRLB

In this section, we consider the misspecified ML (MML) estimation [25, 26] for
a scalar parameter $\theta\in\mathbb{R}$. In order to use the asymptotic
results for the MML estimate $\hat{\theta}$, we are going to assume that the
elements $x_{n}$, $n=0,\ldots,N-1$, of the measurement vector $\mathbf{x}$ are
independent and identically distributed as $x_{n}\sim\bar{f}(x_{n})$ where
$\bar{f}(\cdot)$ denotes the true measurement distribution. The true
distribution for the measurement vector $\mathbf{x}$ is then given as

$\displaystyle\bar{f}(\mathbf{x})=\prod_{n=0}^{N-1}\bar{f}(x_{n}).$ (19)

We assume that MML estimate $\hat{\theta}$ is calculated by maximizing the
assumed likelihood $f(\mathbf{x};\theta)$ given in (15). The relationship of
the suggested variance expression $\widehat{V}_{\hat{\theta}}(\theta)$ to
MCRLB [24, 25] is given in the following proposition.

######  Proposition 2.

  * Assume that

  * A0

The parameter $\theta$ has finite support, i.e.,
$\theta\in[\theta_{\mathrm{min}},\theta_{\mathrm{max}}]$.

  * A1

The first three derivatives of $\ln f(x;\theta)$ with respect to $\theta$
exist for all $\theta$ and are continuous with respect to $\theta$.

  * A2

For every $\theta$, the functions $|\frac{\partial^{i}}{\partial\theta^{i}}\ln
f(x;\theta)|$, $i=0,1,2,3$, are dominated by functions $b_{i}(x)$,
$i=0,1,2,3$, which all have finite variance with respect to the true
measurement distribution $\bar{f}(x)$.

  * A3

The KL divergence $D(\bar{f}(x)||f(x;\theta))$ has a unique minimum with
respect to $\theta$ at
$\theta=\theta_{*}\in(\theta_{\mathrm{min}},\theta_{\mathrm{max}})$.222Note
that the existence of the KL divergence $D(\bar{f}(x)||f(x;\theta))$
necessitates additionally the existence of the $E_{\bar{f}}[\ln\bar{f}(x)]$,
which we implicitly assume for the sake of conceptual simplicity. We may
eliminate the need for the existence of $E_{\bar{f}}[\ln\bar{f}(x)]$ by
stating this assumption differently as in [22, 33, 34].

  * A4

The expectation
$\mathbb{E}_{\bar{f}}\big{[}\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(x;\theta_{*})\big{]}$ is non-zero.

The assumption A0 is sufficient (but not necessary) for the convergence of the
integral in (6). Under the assumptions A1-A3 it can be shown that (See [22,
Theorem 2.1]) the MML estimate $\hat{\theta}\triangleq\arg\max_{\theta}\ln
f(\mathbf{x};\theta)$ is misspecified consistent, i.e.,
$\hat{\theta}\overset{\text{a.s.}}{\rightarrow}\theta_{*}$ as
$N\rightarrow\infty$ . Then, we have

$\displaystyle\widehat{V}_{\mathrm{MML}}(\theta_{*})\rightarrow
V_{\mathrm{MML}}(\theta_{*})$ (20)

as $N\rightarrow\infty$, i.e., the finite support version of the expression
$\widehat{V}_{\mathrm{MML}}(\theta_{*})$ converges to the true large sample
asymptotic variance $V_{\mathrm{MML}}(\theta_{*})$ of the MML estimate
$\hat{\theta}$ as $N\rightarrow\infty$. If the MML estimate is also
asymptotically misspecified efficient, then we have

$\displaystyle\widehat{V}_{\mathrm{MML}}(\theta_{*})\rightarrow\frac{\mathcal{B}(\theta_{*})}{\mathcal{A}^{2}(\theta_{*})}$
(21)

as $N\rightarrow\infty$ where the quantity
$\frac{\mathcal{B}(\theta_{*})}{\mathcal{A}^{2}(\theta_{*})}$ is the MCRLB and

$\mathcal{A}(\theta_{*})\triangleq\,\mathbb{E}_{\bar{f}}\bigg{[}\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(x;\theta_{*})\bigg{]},\quad\mathcal{B}(\theta_{*})\triangleq\,\mathbb{E}_{\bar{f}}\bigg{[}\bigg{(}\frac{\partial}{\partial\theta}\ln
f(x;\theta_{*})\bigg{)}^{2}\bigg{]}.$ (22)

Proof: The proof is given in C. $\blacksquare$

Note that according to the proposition, $\widehat{V}(\theta_{*})$ converges to
the asymptotic variance $V_{\mathrm{MML}}(\theta_{*})$ of the MML estimate
$\hat{\theta}$. When the true measurement distribution $\bar{f}(\cdot)$ admits
the same parameterization as the assumed measurement distribution
$f(\cdot;\theta)$ with the true parameter value $\theta=\bar{\theta}$, i.e.,
$\bar{f}(x)=\bar{f}(x;\bar{\theta})$, then we might predict the MSE
performance of the MML estimator $\hat{\theta}$ as

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{MML}}(\bar{\theta})\triangleq\widehat{V}_{\mathrm{MML}}(\theta_{*})+(\theta_{*}-\bar{\theta})^{2},$
(23)

which would converge to the true MSE of the MML estimator as
$N\rightarrow\infty$ if the assumptions of Proposition 2 are satisfied.

###  4.3 Relationship to ZZB

We consider a Bayesian estimation problem where the parameter $\theta$ is
assigned with the prior distribution $f(\theta)$. The MAP estimate of $\theta$
can then be defined as follows.

$\displaystyle\hat{\theta}\triangleq\arg\max_{\theta}f(\mathbf{x}|\theta)f(\theta)$
(24)

where the likelihood $f(\mathbf{x};\theta)$ is shown with the conditioning
notation as $f(\mathbf{x}|\theta)$ since $\theta$ is now a random variable.
Note that the MAP estimator given above corresponds to an IDE with the
objective function $\mathcal{L}(\mathbf{x};\theta)\triangleq
f(\mathbf{x}|\theta)f(\theta)$. The true MSE of the MAP estimator is given as

$\displaystyle\operatorname{MSE}_{\mathrm{MAP}}\triangleq$
$\displaystyle\int\int(\hat{\theta}-\theta)^{2}f(\mathbf{x}|\theta)\;\mathrm{d}{}\mathbf{x}f(\theta)\;\mathrm{d}{}\theta$
(25a) $\displaystyle=$
$\displaystyle\,\mathbb{E}\big{[}\mathbb{E}\big{[}(\hat{\theta}-\theta)^{2}\big{|}\theta\big{]}\big{]}$
(25b) $\displaystyle=$
$\displaystyle\,\mathbb{E}\big{[}\operatorname{MSE}_{\mathrm{MAP}}(\theta)\big{]},$
(25c)

where the outer expectation in (25b) is with respect to the random variable
$\theta$ and $\operatorname{MSE}_{\mathrm{MAP}}(\theta)$ denotes the true MSE
of the MAP estimator when $\theta$ is given, i.e.,

$\displaystyle\operatorname{MSE}_{\mathrm{MAP}}(\theta)\triangleq\mathbb{E}\big{[}(\hat{\theta}-\theta)^{2}\big{|}\theta\big{]},$
(26)

where the expectation is only with respect to the noisy measurements
$\mathbf{x}$ given $\theta$. Since the problem becomes a non-random parameter
estimation problem when $\theta$ is given, we can predict
$\operatorname{MSE}_{\mathrm{MAP}}(\theta)$ of the MAP estimator using (8) as
follows.

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{MAP}}(\theta)=2\int_{-\infty}^{\infty}\\!|\epsilon|P\Big{(}f(\mathbf{x}|\theta+2\epsilon)f(\theta+2\epsilon)\geq
f(\mathbf{x}|\theta)f(\theta)\Big{|}\theta\Big{)}\;\mathrm{d}{}\epsilon.$ (27)

By substituting the MSE estimate
$\widehat{\operatorname{MSE}}_{\mathrm{MAP}}(\theta)$ in (27) into the place
of $\operatorname{MSE}_{\mathrm{MAP}}(\theta)$ in (25c) we can predict the
overall MSE of the MAP estimate as follows.

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{MAP}}\triangleq\mathbb{E}[\widehat{\operatorname{MSE}}_{\mathrm{MAP}}(\theta)]$
(28a) $\displaystyle\triangleq$
$\displaystyle\,2\int_{-\infty}^{\infty}f(\theta)\int_{-\infty}^{\infty}|\epsilon|P\Big{(}f(\mathbf{x}|\theta+2\epsilon)f(\theta+2\epsilon)\geq
f(\mathbf{x}|\theta)f(\theta)\Big{|}\theta\Big{)}\;\mathrm{d}{}\epsilon\;\mathrm{d}{}\theta$
(28b) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}|\epsilon|f(\theta)P\Big{(}f(\mathbf{x}|\theta+2\epsilon)f(\theta+2\epsilon)\geq
f(\mathbf{x}|\theta)f(\theta)\Big{|}\theta\Big{)}\;\mathrm{d}{}\epsilon\;\mathrm{d}{}\theta$
$\displaystyle+\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}|\epsilon|f(\theta-2\epsilon)P\Big{(}f(\mathbf{x}|\theta-2\epsilon)f(\theta-2\epsilon)\geq
f(\mathbf{x}|\theta)f(\theta)\Big{|}\theta\Big{)}\;\mathrm{d}{}\epsilon\;\mathrm{d}{}\theta$
(28c) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}|\epsilon|f(\theta)P\Big{(}f(\mathbf{x}|\theta+2\epsilon)f(\theta+2\epsilon)\geq
f(\mathbf{x}|\theta)f(\theta)\Big{|}\theta\Big{)}\;\mathrm{d}{}\epsilon\;\mathrm{d}{}\theta$
$\displaystyle+\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}|\epsilon|f(\theta+2\epsilon)P\Big{(}f(\mathbf{x}|\theta)f(\theta)\geq
f(\mathbf{x}|\theta+2\epsilon)f(\theta+2\epsilon)\Big{|}\theta+2\epsilon\Big{)}\;\mathrm{d}{}\epsilon\;\mathrm{d}{}\theta$
(28d) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}|\epsilon|(f(\theta)+f(\theta+2\epsilon))$
$\displaystyle\times\Big{[}\pi_{1}P\Big{(}\pi_{2}f(\mathbf{x}|\theta+2\epsilon)\geq\pi_{1}f(\mathbf{x}|\theta)\Big{|}\theta\Big{)}+\pi_{2}P\Big{(}\pi_{1}f(\mathbf{x}|\theta)\geq\pi_{2}f(\mathbf{x}|\theta+2\epsilon)\Big{|}\theta+2\epsilon\Big{)}\Big{]}\;\mathrm{d}{}\epsilon\;\mathrm{d}{}\theta$
(28e) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}|\epsilon|(f(\theta)+f(\theta+2\epsilon))P_{\min}^{e}(\theta,\theta+2\epsilon)\;\mathrm{d}{}\epsilon\;\mathrm{d}{}\theta$
(28f) $\displaystyle=$
$\displaystyle\,2\int_{0}^{\infty}\int_{-\infty}^{\infty}\epsilon(f(\theta)+f(\theta+2\epsilon))P_{\min}^{e}(\theta,\theta+2\epsilon)\;\mathrm{d}{}\theta\;\mathrm{d}{}\epsilon$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{\infty}\int_{-\infty}^{\infty}\epsilon(f(\theta)+f(\theta+\epsilon))P_{\min}^{e}(\theta,\theta+\epsilon)\;\mathrm{d}{}\theta\;\mathrm{d}{}\epsilon,$
(28g)

where $P_{\min}^{e}(\theta_{1},\theta_{2})$ is the minimum probability of
error for the binary hypothesis testing problem given below.

$\displaystyle\mathcal{H}_{1}:$ $\displaystyle\,\,\mathbf{x}\sim
f(\mathbf{x}|\theta_{1}),$ (29a) $\displaystyle\mathcal{H}_{2}:$
$\displaystyle\,\,\mathbf{x}\sim f(\mathbf{x}|\theta_{2}),$ (29b)

with the prior hypothesis probabilities $P(\mathcal{H}_{1})=\pi_{1}$ and
$P(\mathcal{H}_{2})=\pi_{2}=1-\pi_{1}$ where

$\displaystyle\pi_{1}\triangleq$
$\displaystyle\frac{f(\theta_{1})}{f(\theta_{1})+f(\theta_{2})},$
$\displaystyle\pi_{2}\triangleq$
$\displaystyle\frac{f(\theta_{2})}{f(\theta_{1})+f(\theta_{2})}.$ (30)

The expression (28g) can be seen to be the ZZB (See [7, Eqn. (14)]) without
the so-called valley filling function. As a result
$\widehat{\operatorname{MSE}}_{\mathrm{MAP}}$ calculated using (8) in a
Bayesian framework is equal to the ZZB. Note that this equality is satisfied
irrespective of whether the objective function $f(\mathbf{x}|\theta)f(\theta)$
satisfies the assumptions of Theorem 1 or not. If the objective function
$f(\mathbf{x}|\theta)f(\theta)$, which is actually the joint density
$f(\mathbf{x},\theta)$ of $\mathbf{x}$ and $\theta$, also satisfies the
conditions of Theorem 1, then this would mean that
$\operatorname{MSE}_{\mathrm{MAP}}(\theta)=\widehat{\operatorname{MSE}}_{\mathrm{MAP}}(\theta)$
for all $\theta\in\mathbb{R}$ and hence
$\operatorname{MSE}_{\mathrm{MAP}}=\widehat{\operatorname{MSE}}_{\mathrm{MAP}}=\mathrm{ZZB}$
and hence ZZB would have to be tight, i.e., ZZB would have to be equal to the
true average MSE of the MAP estimate $\hat{\theta}$. As a result, the
conditions of Theorem 1 are also a set of sufficient conditions for ZZB to be
tight.

##  5 Extension to the Case with Nuisance Parameters

Suppose now that we have $J>1$ unknown scalar parameters, i.e.,
$\boldsymbol{\theta}\in\mathbb{R}^{J}$, and we would like to estimate only one
of them while keeping the others as unknown nuisance parameters. Without loss
of generality we assume that we would like to estimate $\theta_{1}$ while
treating the other parameters $\theta_{2},\ldots,\theta_{J}$ as nuisance
parameters. We can express the estimate $\hat{\theta}_{1}$ for $\theta_{1}$ as

$\displaystyle\hat{\theta}_{1}\triangleq[\boldsymbol{\hat{\theta}}]_{1}=$
$\displaystyle\arg\max_{\theta_{1}}\Big{[}\underbrace{\max_{\theta_{\backslash
1}}\mathcal{L}(\mathbf{x};\boldsymbol{\theta})}_{\triangleq\mathcal{L}_{1}(\mathbf{x},\theta_{1})}\Big{]}=\arg\max_{\theta_{1}}\mathcal{L}_{1}(\mathbf{x},\theta_{1}),$
(31)

where $\boldsymbol{\theta}_{\backslash
1}\triangleq\left[\theta_{2}\,\,\theta_{3}\,\cdots\,\theta_{J}\right]^{\mathrm{T}}{}$.
If we assume that the function $\mathcal{L}_{1}(\mathbf{x},\theta_{1})$
defined as
$\mathcal{L}_{1}(\mathbf{x},\theta_{1})\triangleq\max_{\boldsymbol{\theta}_{\backslash
1}}\mathcal{L}(\mathbf{x};\boldsymbol{\theta})$ satisfies the conditions in
Theorem 1, applying the result of Remark 1 to the IDE in (31) would give

$\displaystyle\widehat{\operatorname{MSE}}(\bar{\theta}_{1})=$
$\displaystyle\,2\int_{-\infty}^{\infty}|\epsilon|P\left(\mathcal{L}_{1}(\mathbf{x},\bar{\theta}_{1}+2\epsilon)\geq\mathcal{L}_{1}(\mathbf{x},\bar{\theta}_{1})\right)\;\mathrm{d}{}\epsilon,$
$\displaystyle=$
$\displaystyle\,2\int_{-\infty}^{\infty}|\epsilon|P\Big{(}\max_{\boldsymbol{\theta}_{\backslash
1}}\mathcal{L}(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})\geq\max_{\boldsymbol{\theta}_{\backslash
1}}\mathcal{L}(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\theta}_{\backslash
1})\Big{)}\;\mathrm{d}{}\epsilon.$ (32)

With the selection $\mathcal{L}(\mathbf{x};\boldsymbol{\theta})\triangleq
f(\mathbf{x};\boldsymbol{\theta})\geq 0$, we can obtain the MSE of the ML
estimate $\hat{\theta}_{1}$ of $\theta_{1}$ similarly to Corollary 1 from (32)
as

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta}_{1})=$
$\displaystyle\,2\int_{-\infty}^{\infty}|\epsilon|P\bigg{(}\genfrac{}{}{}{}{\raisebox{2.0pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-2.0pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\theta}_{\backslash 1})$}}\geq
1\bigg{)}\;\mathrm{d}{}\epsilon,$ (33)

connecting the MSE of the ML estimator to the error probability of a
generalized likelihood ratio test (GLRT) (instead of a likelihood ratio test)
in the presence of nuisance parameters [35].

In Section 6.2 below, we are going to investigate the expression (33) further
on the specific case of the parametric mean model and make approximations to
facilitate its calculation, which are later extended to the general case in a
remark.

##  6 Application to ML Estimation with the Parametric Mean Model with
Gaussian Noise

We consider ML estimator with the measurement model given as

$\displaystyle\mathbf{x}=\mathbf{m}(\boldsymbol{\bar{\theta}})+\mathbf{v},$
(34)

where
$\mathbf{v}\sim\mathcal{CN}(\mathbf{v};\mathbf{0},\sigma^{2}\mathbf{I}_{N})$
represents the measurement noise and the manifold function
$\mathbf{m}:\mathbb{R}^{J}\rightarrow\mathbb{C}^{N}$ is, in general, a
complex-valued function of the unknown parameter vector
$\boldsymbol{\theta}\in\mathbb{R}^{J}$. The measurement model in (34) is
widely used in signal processing applications. For example, a linear manifold
function
$\mathbf{m}(\boldsymbol{\bar{\theta}})=\mathbf{H}\boldsymbol{\bar{\theta}}$
may represent a multi-input multi-output (MIMO) communication system; a non-
linear manifold function may represent the array response in the direction of
arrival estimation problems [1].

The likelihood function for an arbitrary $\boldsymbol{\theta}$ is given as

$\displaystyle
f(\mathbf{x};\boldsymbol{\theta})=\mathcal{CN}(\mathbf{x};\mathbf{m}(\boldsymbol{\theta}),\sigma^{2}\mathbf{I}_{N}).$
(35)

We investigate the cases of a scalar parameter with and without nuisance
parameters in different subsections below. In order to calculate the predicted
MSE values we will need the following log-likelihood ratio expression.

$\displaystyle\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\boldsymbol{\theta})$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\boldsymbol{\bar{\theta}})$}}=$
$\displaystyle\frac{1}{\sigma^{2}}\big{(}2\Re\\{\mathbf{\tilde{m}}^{\mathrm{H}}{}(\boldsymbol{\theta};\boldsymbol{\bar{\theta}})(\mathbf{x}-\mathbf{m}(\bar{\theta}))\\}-\|\mathbf{\tilde{m}}(\boldsymbol{\theta};\boldsymbol{\bar{\theta}})\|^{2}\big{)},$
(36)

where
$\mathbf{\tilde{m}}(\boldsymbol{\theta}_{1};\boldsymbol{\theta}_{2})\triangleq\mathbf{m}(\boldsymbol{\theta}_{1})-\mathbf{m}(\boldsymbol{\theta}_{2}).$

###  6.1 Case of a Scalar Parameter with No Nuisance Parameters

Suppose now that we have a scalar parameter $\theta$ with the true value
$\bar{\theta}$ ($J=1$). Note that this case can also be interpreted to be the
case when we have multiple parameters
$\boldsymbol{\theta}=[\theta_{1}\,\,\,\boldsymbol{\theta}_{\backslash 1}]$ and
the true values $\boldsymbol{\bar{\theta}}_{\backslash 1}$ of the nuisance
parameters $\boldsymbol{\theta}_{\backslash 1}$ are perfectly known. We can
evaluate the probability in (13) as

$\displaystyle
P\left(\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
1\right)$
$\displaystyle=P\left(\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
0\right)=P\left(2\Re\\{\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}+2\epsilon;\bar{\theta})\mathbf{v}\\}\geq\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon;\bar{\theta})\|^{2}\right)$
(37a)
$\displaystyle=\mathcal{N}_{\mathrm{\mathrm{ccdf}}}\left(\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon;\bar{\theta})\|^{2};0,2\sigma^{2}\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon;\bar{\theta})\|^{2}\right)$
(37b)
$\displaystyle=\mathcal{N}_{\mathrm{\mathrm{ccdf}}}\left(\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon;\bar{\theta})\|;0,2\sigma^{2}\right),$
(37c)

under the assumption that
$\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon;\bar{\theta})\|\neq 0$, where
$\mathcal{N}_{\mathrm{\mathrm{ccdf}}}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma})$
denotes the complementary cumulative distribution function
($\mathrm{\mathrm{ccdf}}$) of a real Gaussian random vector with mean
$\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$ evaluated at
$\mathbf{x}$. Assuming that
$\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon;\bar{\theta})\|\neq 0$ for almost
all $\epsilon\in\mathbb{R}$, we can substitute this probability expression
into (13) to get

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})=\,2\int_{-\infty}^{\infty}|\epsilon|\mathcal{N}_{\mathrm{ccdf}}\left(\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon,\bar{\theta})\|;0,2\sigma^{2}\right)\;\mathrm{d}{}\epsilon.$
(38)

######  Remark 2.

If both the function $\mathbf{m}(\cdot)$ and the measurement noise
$\mathbf{v}$ are real-valued, i.e., if we have
$\mathbf{m}:\mathbb{R}\rightarrow\mathbb{R}^{N}$ and
$\mathbf{v}\sim\mathcal{N}(\mathbf{v};\mathbf{0},\sigma^{2}\mathbf{I}_{N})$,
then, instead of (38), one needs to use

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})=\,2\int_{-\infty}^{\infty}|\epsilon|\mathcal{N}_{\mathrm{ccdf}}\left(\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon,\bar{\theta})\|;0,4\sigma^{2}\right)\;\mathrm{d}{}\epsilon.~{}$
(39)

Using this expression amounts to replacing the variance $\sigma^{2}$ in (38)
with $2\sigma^{2}$. $\blacksquare$

Note that the likelihood (35) does not satisfy the conditions of Theorem 1 and
its corollary in general except for some trivial cases, e.g., the case of
linear or affine manifold function $\mathbf{m}(\theta)$. As a result, the
predicted MSE expressions in (38) and (39) are expected to be only an
approximate estimate of the true MSE of the ML estimator. Furthermore, a
closed form solution rarely exists for the integrals in (38) and (39).
Therefore, numerical integration methods have to be used as shown in Example 1
below.

The relations in (38) and (39) provide some insight on the suggested MSE
expression. As $\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon,\bar{\theta})\|$,
which is the norm of the difference between manifold vectors
$\mathbf{m}(\bar{\theta}+2\epsilon)$ and $\mathbf{m}(\bar{\theta})$, gets
larger, it should be easier to accurately estimate $\theta$ and we get a
smaller predicted MSE value (since the function
$\mathcal{N}_{\mathrm{ccdf}}(\cdot)$ monotonically decreases as its argument
gets larger). Also, it is interesting to see that, for the simplest case
$\mathbf{{m}(\theta)}\triangleq\theta$, the expression for the predicted MSE
in (38) simplifies to;

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})$
$\displaystyle=\,2\int_{-\infty}^{\infty}|\epsilon|\mathcal{N}_{\mathrm{ccdf}}\left(|2\epsilon|;0,2\sigma^{2}\right)\;\mathrm{d}{}\epsilon=\frac{\sigma^{2}}{2},$
(40)

which can be obtained by integration by parts. This result is expected, as the
corresponding ML estimator is $\hat{\theta}=\Re\\{x\\}$, hence, the true MSE
must be equal to half of the noise variance. We have an exact result since the
objective function is a Gaussian likelihood satisfying the conditions of
Theorem 1. We finally consider the following example in order to illustrate
the practical simplicity of the expression (38).

######  Example 1 (Frequency estimation using ML).

Consider the following signal model.

$\displaystyle x_{n}=Ae^{j\bar{\omega}n}+v_{n},\quad n=0,\ldots,N-1,$ (41)

where $A\in\mathbb{C}$ is the known complex amplitude;
$\bar{\omega}\in[-\pi,\pi]$ is the unknown true frequency to be estimated
using the ML estimator; $v_{n}\sim\mathcal{CN}(v_{n};0,\sigma^{2})$,
$n=0,\ldots,N-1$, is the white measurement noise. MSE of the ML estimator
based on the measurements $x_{n}$, $n=0,\ldots,N-1$, can be calculated with
the Matlab function given in Figure 2, which involves only three lines of
code. A sample run can be made using the command _ML_frequency(pi/2,1,1,16)
for the true frequency value $\bar{\omega}=\pi/2$, amplitude $A=1$, noise
variance $\sigma^{2}=1$ and number of samples $N=16$ gives
$\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\frac{\pi}{2})=$17e-4
$\mathrm{rad}^{2}$. Note that the function culateMSEhat(.) in Figure 2 can be
used for predicting the MSE performance of the ML estimator for any
measurement model of type (34) for a scalar parameter
$\theta\in[\theta_{\mathrm{min}},\theta_{\mathrm{max}}]$.

⬇

1ction MSEhatML = MSE_ML_frequency(tw,A,sigma2,N)

2% Calculate MSEhatML for ML frequency estimate

3% tw: true value of the frequency (1x1) (rad/sec) (-pi < tw < pi)

4% A: known complex amplitude (1x1)

5% sigma2: known complex normal noise variance (1x1)

6% N: number of samples (1x1)

7m = @(w) A*exp (1i*(0:N-1)'*w);

8MSEhatML = calculateMSEhat(m,tw,-pi,pi,sigma2);

9

10function MSEhatML = calculateMSEhat(m,tt,tmin,tmax,sigma2)

11% Calculates the MSEhat for the parameterized mean model

12% m: function handle: m(.) takes a 1xL array of parameter values [theta1
theta2 ... thetaL]

13% : and it returns the Nm x L matrix [m(theta1) m(theta2) ... m(thetaL)]

14% tt: true theta value (1x1)

15% tmin: minimum value of theta (1x1)

16% tmax: maximum value of theta (1x1)

17% sigma2: complex normal noise variance (1x1)

18MSEhatML =
2*integral(@(e)abs(e).*(1-normcdf(vecnorm(m(tt+2*e)-m(tt),2,1)/sqrt(2*sigma2))),...

19 (tmin-tt)/2,(tmax-tt)/2);

Figure 2: A Matlab code (R2021b) for predicting the MSE of the ML estimator
for the frequency estimation problem.

###  6.2 Case of a Scalar Parameter with Nuisance Parameters

When some nuisance parameters exist, we consider the case in Section 5 and use
the MSE expression in (33). Unfortunately it is analytically difficult to
calculate the maxima and the probabilities in the integrands on the right hand
side of (33) exactly. In the following, we are going to make some
approximations to facilitate the calculation. Similar approximations can also
be made for the more general case in (32) (See Remark 3 below).

$\displaystyle\genfrac{}{}{}{}{\raisebox{1.5pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-1.5pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\theta}_{\backslash
1})$}}\approx$
$\displaystyle\,\genfrac{}{}{}{}{\raisebox{1.5pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})$}}=\max_{\boldsymbol{\theta}_{\backslash
1}}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})$}}\approx\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})$}},$ (42)

where the first approximation in (42) is made by assuming that the maximum in
the denominator of the left hand side is achieved approximately at the true
values of the nuisance parameters, i.e., at $\boldsymbol{\theta}_{\backslash
1}=\boldsymbol{\bar{\theta}}_{\backslash 1}$, which is reasonable under
asymptotic conditions. The set $\boldsymbol{\Theta}_{\backslash
1}\triangleq\\{\boldsymbol{\theta}_{\backslash
1}^{1},\boldsymbol{\theta}_{\backslash
1}^{2},\ldots,\boldsymbol{\theta}_{\backslash 1}^{N_{\theta}}\\}$ appearing in
(42) is a set of grid points including the true value
$\boldsymbol{\bar{\theta}}_{\backslash 1}$ of $\boldsymbol{\theta}_{\backslash
1}$. Using these approximations, we can approximate the probability in (33) as

$\displaystyle P\bigg{(}$
$\displaystyle\genfrac{}{}{}{}{\raisebox{1.5pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-1.5pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\theta}_{\backslash 1})$}}\geq
1\bigg{)}\approx\,P\bigg{(}\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})$}}\geq 1\bigg{)},$ (43a) $\displaystyle=$
$\displaystyle\,P\bigg{(}\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})$}}\geq 0\bigg{)}$ (43b) $\displaystyle=$
$\displaystyle\,P\bigg{(}\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\Big{[}2\Re\\{\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\mathbf{v}\\}-\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})\|^{2}\Big{]}\geq
0\bigg{)}$ (43c) $\displaystyle=$
$\displaystyle\,1-P\bigg{(}\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\Big{[}2\Re\\{\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\mathbf{v}\\}-\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})\|^{2}\Big{]}\leq
0\bigg{)}$ (43d)

where $\mathbf{\tilde{m}}(\theta_{1}^{1},\boldsymbol{\theta}_{\backslash
1}^{1};\theta_{1}^{2},\boldsymbol{\theta}_{\backslash
1}^{2})\triangleq\mathbf{m}(\theta_{1}^{1},\boldsymbol{\theta}_{\backslash
1}^{1})-\mathbf{m}(\theta_{1}^{2},\boldsymbol{\theta}_{\backslash 1}^{2}).$

Let us now define the matrix
$\mathbf{\widetilde{M}}_{\epsilon}\in\mathbb{R}^{N\times N_{\theta}}$ and the
vector $\boldsymbol{\tilde{\mu}}_{\epsilon}\in\mathbb{R}^{N_{\theta}}$ as

$\mathbf{\widetilde{M}}_{\epsilon}\triangleq\left[\begin{array}[]{c}\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1}^{1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})\\\ \vdots\\\
\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1}^{N_{\theta}};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\end{array}\right]^{\mathrm{H}}{},\qquad\boldsymbol{\tilde{\mu}}_{\epsilon}\triangleq\left[\begin{array}[]{c}\big{\|}\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1}^{1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\big{\|}^{2}\\\ \vdots\\\
\big{\|}\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1}^{N_{\theta}};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\big{\|}^{2}\end{array}\right].$ (44)

We can now write (43d) as

$P\bigg{(}\genfrac{}{}{}{}{\raisebox{1.5pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-1.5pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\theta}_{\backslash 1})$}}\geq
1\bigg{)}\approx\,1-P\Big{(}2\Re\\{\mathbf{\widetilde{M}}_{\epsilon}^{\mathrm{H}}{}\mathbf{v}\\}\leq\boldsymbol{\tilde{\mu}}_{\epsilon}\Big{|}\boldsymbol{\bar{\theta}}\Big{)}=\,\mathcal{N}_{\mathrm{ccdf}}\big{(}\boldsymbol{\tilde{\mu}}_{\epsilon};\mathbf{0},2\sigma^{2}\Re\big{\\{}\mathbf{\widetilde{M}}_{\epsilon}^{\mathrm{H}}{}\mathbf{\widetilde{M}}_{\epsilon}\big{\\}}\big{)},$
(45)

under the assumption that
$\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1}^{i};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})\|\neq 0$ for
$i=1,\ldots,N_{\theta}$. Note that the inequalities between vector quantities
above should be interpreted in an elementwise manner. The probability given in
(45) is the generalization of the single parameter probability in (37b) to the
case of (presence of) nuisance parameters. In fact, when we select the set
$\boldsymbol{\Theta}_{\backslash 1}$ as $\boldsymbol{\Theta}_{\backslash
1}=\\{\boldsymbol{\bar{\theta}}_{\backslash 1}\\}$, i.e., when we have a grid
composed of only the true nuisance parameter
$\boldsymbol{\bar{\theta}}_{\backslash 1}$, the probability (45) reduces to
the probability (37b). Moreover since the set $\boldsymbol{\Theta}_{\backslash
1}$ contains the true value $\boldsymbol{\bar{\theta}}_{\backslash 1}$, the
probability (45) is always larger than or equal to the probability (37b).
Substituting the result (45) into (33) we get

$\displaystyle\widehat{\operatorname{MSE}}(\bar{\theta}_{1})=\,2\int_{-\infty}^{\infty}|\epsilon|\mathcal{N}_{\mathrm{ccdf}}\big{(}\boldsymbol{\tilde{\mu}}_{\epsilon};\mathbf{0},2\sigma^{2}\Re\big{\\{}\mathbf{\widetilde{M}}_{\epsilon}^{\mathrm{H}}{}\mathbf{\widetilde{M}}_{\epsilon}\big{\\}}\big{)}\;\mathrm{d}{}\epsilon.$
(46)

Note that since the probability (45) is always larger than or equal to the
probability (37b), the predicted MSE in (46) is always larger than or equal to
the single parameter predicted MSE in (38).

Although we ended up with an analytical expression for the predicted MSE in
the nuisance parameter case, unfortunately, the calculation of the predicted
MSE in (46) involves the numerical calculation of the $N_{\theta}$-variate
normal (c)cdf which can be carried out for only small values of the number of
grid points $N_{\theta}$. Furthermore, the covariance matrix
$\Re\big{\\{}\mathbf{\widetilde{M}}_{\epsilon}^{\mathrm{H}}{}\mathbf{\widetilde{M}}_{\epsilon}\big{\\}}$
might be ill-conditioned or singular which makes the calculation of the
probability even more difficult. As a result, the calculation of the predicted
MSE in (46) would be computationally infeasible for large grid sizes
$N_{\theta}$. To avoid this calculation we might follow an alternative
approach by approximating the right hand side of (43c) as

$\displaystyle P\bigg{(}$
$\displaystyle\genfrac{}{}{}{}{\raisebox{1.5pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x}|\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})$}}{\raisebox{-1.5pt}{$\max_{\boldsymbol{\theta}_{\backslash
1}}f(\mathbf{x}|\bar{\theta}_{1},\boldsymbol{\theta}_{\backslash 1})$}}\geq
1\bigg{)}$
$\displaystyle\approx\,P\bigg{(}\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\Big{[}2\Re\\{\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\mathbf{v}\\}-\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})\|^{2}\Big{]}\geq
0\bigg{)}$ (47a) $\displaystyle\approx\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}P\bigg{(}\Big{[}2\Re\\{\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\mathbf{v}\\}-\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})\|^{2}\Big{]}\geq
0\bigg{)}$ (47b) $\displaystyle=\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}P\Big{(}2\Re\\{\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\mathbf{v}\\}\geq\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})\|^{2}\Big{)}$
(47c) $\displaystyle=\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\mathcal{N}_{\mathrm{ccdf}}\Big{(}\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\|^{2};0,2\sigma^{2}\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})\|^{2}\Big{)}$
(47d) $\displaystyle=\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\mathcal{N}_{\mathrm{ccdf}}\big{(}\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\|;0,2\sigma^{2}\big{)}$ (47e)
$\displaystyle=\mathcal{N}_{\mathrm{ccdf}}\Big{(}\min_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\|;0,2\sigma^{2}\Big{)},$ (47f)

under the assumption that $\min_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})\|\neq 0$. The
approximation sign in (47a) represents the approximations made until reaching
(43c). The approximation sign in (47b) can be replaced with a greater than
equal to sign, i.e., the right hand side of it is a lower bound for the left
hand side. Substituting (47f) into (33) gives the predicted MSE expression
shown below:

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta}_{1})=$
$\displaystyle\,2\int_{-\infty}^{\infty}|\epsilon|\mathcal{N}_{\mathrm{ccdf}}\Big{(}\min_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\|;0,2\sigma^{2}\Big{)}\;\mathrm{d}{}\epsilon.$ (48)

The predicted MSE in (48) is always smaller than or equal to the
computationally prohibitive predicted MSE in (46) due to the approximation
made in (47b), however, it requires the calculation of the ccdf of only a
univariate normal random variable. Note that the MSE in (48) is still always
larger than or equal to the single parameter MSE in (38) since the parameter
grid $\boldsymbol{\Theta}_{\backslash 1}$ contains the true value
$\boldsymbol{\bar{\theta}}_{\backslash 1}$ of $\boldsymbol{\theta}_{\backslash
1}$. This is because of the fact that

$\displaystyle\min_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash 1}}$
$\displaystyle\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\|\leq\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\bar{\theta}}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\|=\|\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon;\bar{\theta}_{1})\|$
(49)

and that the function $\mathcal{N}_{\mathrm{ccdf}}(\cdot,0,2\sigma^{2})$
monotonically increases as its argument gets smaller.

The intuitive meaning of the MSE expression (48) can be explained as follows.
When the nuisance parameters $\boldsymbol{\theta}_{\backslash 1}$ are known,
i.e., we have the case of a single parameter in Section 6.1, the MSE is seen
to be dependent on the distance between the mean vector
$\mathbf{m}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\bar{\theta}}_{\backslash
1})$ and the true mean vector
$\mathbf{m}(\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})$, which
was shown as (the magnitude of) the vector
$\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon;\bar{\theta}_{1})\triangleq\mathbf{\tilde{m}}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\bar{\theta}}_{\backslash
1};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})$ in (38). On the
other hand, when the the nuisance parameters $\boldsymbol{\theta}_{\backslash
1}$ are not known, the predicted MSE is dependent on minimum distance between
the mean vectors
$\mathbf{m}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash 1})$,
where $\boldsymbol{\theta}_{\backslash 1}$ takes values in a grid containing
the true nuisance parameter value $\boldsymbol{\bar{\theta}}_{\backslash 1}$,
and the fixed true mean vector
$\mathbf{m}(\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})$. Hence
if the vector
$\mathbf{m}(\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash 1})$ is
similar to the fixed vector
$\mathbf{m}(\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash 1})$ for
some values of the nuisance parameter $\boldsymbol{\theta}_{\backslash 1}$ in
the grid, the resulting predicted MSE would get larger.

######  Remark 3.

The approximations made in this section on (33) can be applied to the general
case (32) as follows.

$\displaystyle P\Big{(}\max_{\boldsymbol{\theta}_{\backslash
1}}\mathcal{L}(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})\geq\max_{\boldsymbol{\theta}_{\backslash
1}}\mathcal{L}(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\theta}_{\backslash
1})\Big{)}\approx$
$\displaystyle\,P\Big{(}\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}\mathcal{L}(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})\geq\mathcal{L}(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\Big{)},$ (50a) $\displaystyle\approx$
$\displaystyle\max_{\boldsymbol{\theta}_{\backslash
1}\in\boldsymbol{\Theta}_{\backslash
1}}P\big{(}\mathcal{L}(\mathbf{x};\bar{\theta}_{1}+2\epsilon,\boldsymbol{\theta}_{\backslash
1})\geq\mathcal{L}(\mathbf{x};\bar{\theta}_{1},\boldsymbol{\bar{\theta}}_{\backslash
1})\big{)}.$ (50b)

###  6.3 Application to ML Estimation under Model Mismatch

In this section we consider the problem of ML estimation with the parametric
mean model under model mismatch, also known as misspecified ML (MML)
estimation in the literature [26, 25]. For the sake of simplicity we consider
only the scalar parameter case, i.e., $\theta\in\mathbb{R}$. The measurements
$\mathbf{x}$ are modeled as
$\mathbf{x}=\mathbf{\bar{m}}(\bar{\theta})+\mathbf{v}$, where
$\mathbf{\bar{m}}(\cdot)$ denotes the true mean function and
$\mathbf{v}\sim\mathcal{CN}(\mathbf{v};\mathbf{0},\bar{\sigma}^{2}\mathbf{I}_{N})$
represents the measurement noise with the true variance $\bar{\sigma}^{2}$.
This model corresponds to the true likelihood
$\bar{f}(\mathbf{x};\bar{\theta})\triangleq\mathcal{CN}(\mathbf{x};\mathbf{\bar{m}}(\bar{\theta}),\bar{\sigma}^{2}\mathbf{I}_{N}).$
We are interested in the MSE of the mismatched ML estimator $\hat{\theta}$ of
$\theta$ given as

$\displaystyle\hat{\theta}\triangleq\arg\max_{\theta}f(\mathbf{x};\theta),$
(51)

where the objective function is the assumed likelihood $f(\mathbf{x};\theta)$
given as

$\displaystyle
f(\mathbf{x};\theta)\triangleq\mathcal{CN}(\mathbf{x};\mathbf{m}(\theta),\sigma^{2}\mathbf{I}).$
(52)

For predicting the performance of the MML estimator given above, we can use
the MSE expression of Remark 1 by setting
$\mathcal{L}(\mathbf{x};\theta)\triangleq f(\mathbf{x};\theta)$ and
calculating the probability in the integrand of (8) with respect to the true
measurement distribution $\bar{f}(\cdot;\bar{\theta})$. The log-likelihood
ratio
$\ln\frac{f(\mathbf{x};\bar{\theta}+2\epsilon)}{f(\mathbf{x};\bar{\theta})}$
in this case is given as

$\displaystyle\ln$
$\displaystyle\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}=\frac{1}{\sigma^{2}}\big{(}2\Re\big{\\{}\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}+2\epsilon;\bar{\theta})(\mathbf{x}-\mathbf{\bar{m}}(\bar{\theta}))\big{\\}}-\|\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon;\bar{\theta})\|^{2}+2\Re\big{\\{}\mathbf{\tilde{m}}^{\mathrm{H}}{}(\bar{\theta}+2\epsilon;\bar{\theta})\boldsymbol{\mu}(\bar{\theta})\big{\\}}\big{)}$
(53)

where $\mathbf{\tilde{m}}(\cdot,\cdot)$ was defined in (6) and

$\displaystyle\boldsymbol{\mu}(\theta)\triangleq\mathbf{\bar{m}}(\theta)-\mathbf{m}(\theta).$
(54)

We can now calculate the probability of the event
$\ln({f(\mathbf{x};\bar{\theta}+2\epsilon)}/{f(\mathbf{x};\bar{\theta})})\geq
0$ with respect to the true measurement distribution
$\bar{f}(\cdot;\bar{\theta})$ as

$\displaystyle P$
$\displaystyle\bigg{(}\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x}|\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x}|\bar{\theta})$}}\geq
0\bigg{)}=\mathcal{N}_{\mathrm{ccdf}}\bigg{(}\|\mathbf{\tilde{m}}(\cdot)\|-2\Re\bigg{\\{}\frac{\mathbf{\tilde{m}}^{\mathrm{H}}{}(\cdot)}{\|\mathbf{\tilde{m}}(\cdot)\|}\boldsymbol{\mu}(\bar{\theta})\bigg{\\}};0,2\bar{\sigma}^{2}\bigg{)}$
(55)

under the assumption that $\|\mathbf{\tilde{m}}(\cdot)\|\neq 0$, where we
dropped the arguments of the function
$\mathbf{\tilde{m}}(\bar{\theta}+2\epsilon;\bar{\theta})$ for brevity.
Substituting this expression into the integrand of (8) we get the following
predicted MSE for the MML estimate.

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{MML}}$
$\displaystyle(\bar{\theta})=\,2\int_{-\infty}^{\infty}|\epsilon|\mathcal{N}_{\mathrm{ccdf}}\bigg{(}\|\mathbf{\tilde{m}}(\cdot)\|-2\Re\bigg{\\{}\frac{\mathbf{\tilde{m}}^{\mathrm{H}}{}(\cdot)}{\|\mathbf{\tilde{m}}(\cdot)\|}\boldsymbol{\mu}(\bar{\theta})\bigg{\\}};0,2\bar{\sigma}^{2}\bigg{)}\;\mathrm{d}{}\epsilon.$
(56)

##  7 Numerical Results

In this section, we examine the performance of the proposed MSE expression on
four different direction of arrival (DOA) estimation problems. The first two
problems study the conventional and misspecified ML estimation respectively.
In the third one, we investigate the performance of an IDE whose objective
function is not the likelihood function, but a function derived from the
manifold characteristics. The fourth problem investigates Bayesian DOA
estimation. The implementation details of the numerical experiments are given
in D.

###  7.1 DOA Estimation (No Model Mismatch)

Consider the DOA estimation problem with an $N$-element sensor array with the
following array manifold.

$\mathbf{a}_{\boldsymbol{\psi}}=\,[a_{1},a_{2},\ldots,a_{N}]^{\mathrm{T}}{},\quad
a_{n}=\exp\left(j\frac{2\pi}{\lambda}\mathbf{p}_{n}^{\mathrm{T}}{}\mathbf{u}_{\boldsymbol{\psi}}\right),\quad\mathbf{u}_{\boldsymbol{\psi}}=\left[\begin{array}[]{c}\cos(\phi)\sin(\theta)\\\
\sin(\phi)\sin(\theta)\\\
\cos(\theta)\end{array}\right],\quad\mathbf{p}_{n}=\left[\begin{array}[]{c}p_{n}^{\mathsf{x}}\\\
p_{n}^{\mathsf{y}}\\\ p_{n}^{\mathsf{z}}\end{array}\right],$ (57)

where, $\boldsymbol{\psi}\triangleq[\phi,\,\,\theta]^{\mathrm{T}}{}$ denotes
the unknown DOA vector composed of azimuth $\phi\in[0,\,2\pi)$ rads (measured
from the x-axis in counter-clockwise direction) and elevation
$\theta\in[0,\,\pi)$ rads (measured from the z-axis). $\mathbf{p}_{n}$ is the
position vector of the $n$th sensor containing the $\mathsf{x}$, $\mathsf{y}$
and $\mathsf{z}$-coordinates; $N=11$ is the number of sensors; $\lambda$
denotes the wavelength.

Table 1: Sensor positions for the array in Figure 3a Sensor-ID | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11
---|---|---|---|---|---|---|---|---|---|---|---
x, [$\lambda$] | 1.6667 | 1.1785 | 0 | -1.1785 | -1.6667 | -1.1785 | 0 | 1.1785 | 0 | 0 | 0
y, [$\lambda$] | 0 | 1.1785 | 1.6667 | 1.1785 | 0 | -1.1785 | -1.6667 | -1.1785 | 0 | 0 | 0
z, [$\lambda$] | 0 | 1.1785 | 1.6667 | 1.1785 | 0 | -1.1785 | -1.6667 | -1.1785 | 0 | 1.6667 | -1.6667

(a) Array configuration.

(b) Array beampattern at true DOA, $\phi=25^{\circ}$, $\theta=60^{\circ}$.

Figure 3: Array configuration and array beampattern at true DOA.

(a) Azimuth estimation at the true DOA: $\bar{\phi}=25^{\circ}$,
$\bar{\theta}=60^{\circ}$.

(b) Elevation estimation at the true DOA: $\bar{\phi}=25^{\circ}$,
$\bar{\theta}=60^{\circ}$.

Figure 4: Azimuth and elevation estimation performance curves for the array
configuration in Figure 3a. Blue and orange colors correspond to the cases
without and with the nuisance parameter, respectively.

The array, whose sensor positions are given in Table 1, is illustrated in
Figure 3a. The sensor measurement vector $\mathbf{x}\in\mathbb{C}^{N}$ under
additive noise is modeled as

$\mathbf{x}=\beta\mathbf{a}_{\boldsymbol{\bar{\psi}}}+\mathbf{v},$ (58)

where
$\mathbf{v}\sim\mathcal{CN}(\mathbf{v};\mathbf{0},\sigma^{2}\mathbf{I}_{N})$;
$\beta\in\mathbb{R}$ ($\beta$ is taken as a real-valued scalar with no loss of
generality due to the circular symmetry of the complex Gaussian noise) and
$\boldsymbol{\bar{\psi}}\triangleq[\bar{\phi},\,\,\bar{\theta}]^{\mathrm{T}}{}$
denotes the true value of the angle vector $\boldsymbol{\psi}$. The true
target angular positions are $\bar{\phi}=25^{\circ}$ and
$\bar{\theta}=60^{\circ}$. The beampattern of the array, obtained using the
conventional, i.e., Bartlett, beamformer with coefficients steered to the true
DOA, for this angular position is shown in Figure 3b. The beampattern contains
sidelobes as high as $-2$ dB, with a response normalized to $0$ dB at the true
DOA. Consequently the array is prone to gross errors. With these definitions,
the ML estimator involves the following optimization problem:

$\displaystyle\boldsymbol{\hat{\psi}}=$
$\displaystyle\arg\max_{\boldsymbol{\psi}}\Re\\{\mathbf{x}^{\mathrm{H}}{}\mathbf{a}_{\boldsymbol{\psi}}\\}.$
(59)

Since the signal model is a parametric mean model with the mean function
$\mathbf{m}(\boldsymbol{\psi})=\beta\boldsymbol{a_{\psi}}$, the finite support
versions of the expressions (38) and (48) can be utilized for MSE prediction.
We consider three different cases: (i) Azimuth $\phi$ is unknown, but
elevation $\theta=\bar{\theta}$ is known; (ii) Elevation $\theta$ is unknown,
but azimuth $\phi=\bar{\phi}$ is known; (iii) Both azimuth $\phi$ and
elevation $\theta$ are unknown. The results of $10^{5}$ Monte Carlo
simulations are given in Figure 4a and Figure 4b for azimuth and elevation
estimates, respectively. For comparison purposes the corresponding CRLBs (see
[36] for the analytical expressions), Barankin bounds (BBs) [14] with single
test point optimized over a grid, Fessler’s method [27], So et al.’s method
[28], and method of interval errors (MIE) [29] are also illustrated. As seen
from Figure 4, the proposed method is able to predict the threshold
$\operatorname{SNR}$ below which the ML estimator starts following the CRLB
and tracks the CRLB in the asymptotic region as expected. BB, on the other
hand, converges to CRLB at a much smaller SNR value than the ML estimator. MIE
closely follows the ML estimator in the threshold region. This is essentially
due to the problem specific selection of the intervals and accurate gross
error probability calculation. Note that MIE does not have any assumptions on
the objective function, such as symmetry or unimodality, leading to a better
tracking of ML estimator performance especially in the threshold region.
Taylor expansion based methods of Fessler and So et al. follow the CRLB values
in all regions of operation and they are unable to take into account the gross
errors the ML estimator makes below the threshold SNR.

###  7.2 DOA Estimation (Model Mismatch)

In this section we consider the misspecified ML estimation problem examined in
Section 6.3 on the parameterized mean model. For this purpose we consider the
near field azimuth estimation problem with known elevation angle, in which the
estimator uses the plane wave propagation assumption (far field assumption)
rather than the true propagation model which is the spherical spreading.

A uniform circular array of radius $5\lambda/3$ with 12 elements is used. The
signal of interest emanates from a target at a range of $5\lambda$, which is
closer than the far-field limit $2(10\lambda/3)^{2}/\lambda=200\lambda/9$
[37]. The array configuration and the target position are illustrated in
Figure 6.

The true signal model is given as
$\mathbf{\bar{m}}(\bar{\phi})=\,[\bar{a}_{1},\bar{a}_{2},\ldots,\bar{a}_{N}]^{\mathrm{T}}{},\bar{a}_{n}=\exp\big{(}-j\frac{2\pi}{\lambda}d_{n}(\bar{\phi})\big{)},$
where $d_{n}(\phi)=\,\|\mathbf{p}_{n}-r\mathbf{u}_{\phi}\|$,
$\mathbf{u}_{\phi}=\,[\cos(\phi),\,\sin(\phi)]^{\mathrm{T}}{}$,
$\mathbf{p}_{n}=\,[p_{n}^{\mathsf{x}},\,p_{n}^{\mathsf{y}}]^{\mathrm{T}}{}$
and $r$ is the range of the target from the array center as illustrated in
Figure 6. The assumed model by the estimator is the plane wave model, given as
$\mathbf{m}(\phi)=\,[a_{1},a_{2},\ldots,a_{N}]^{\mathrm{T}}{}$,
$a_{n}=\exp\big{(}j\frac{2\pi}{\lambda}\mathbf{p}_{n}^{\mathrm{T}}{}\mathbf{u}_{\phi}\big{)}.$
There is no misspecification in the noise variance, i.e.,
$\sigma^{2}=\bar{\sigma}^{2}$.

The MSE values for this experiment with 10,000 Monte Carlo runs are given in
Figure 6 along with the corresponding MCRLBs [24, 25], BBs [14], and the
results for Fessler’s [27] and So et al.’s [28] methods. MCRLB reduces to the
following expression for this specific problem.

$\displaystyle\operatorname{MCRLB}(\bar{\phi})=\boldsymbol{\operatorname{C}}^{-1}_{D}(\bar{\phi})\mathcal{I}_{D}(\bar{\phi})\boldsymbol{\operatorname{C}}^{-1}_{D}(\bar{\phi}),$
(60)

where

$\displaystyle\mathcal{I}_{D}(\phi)$
$\displaystyle=\frac{2}{\sigma^{2}}\bigg{\|}\frac{\partial\mathbf{m}(\phi)}{\partial\phi}\bigg{\|}^{2},\quad\boldsymbol{\operatorname{C}}_{D}(\theta)=-\mathcal{I}_{D}(\phi)+2\Re\bigg{\\{}\bigg{[}\frac{\partial^{2}\mathbf{m}(\phi)}{\partial\theta^{2}}\bigg{]}^{\mathrm{H}}{}\boldsymbol{\mu}(\phi)\bigg{\\}},$
(61)

and $\boldsymbol{\mu}(\cdot)$ was defined in (54). BB [14] with a single test
point optimized over a grid, which is also the HCRB [12, 15] can be expressed
as follows.

$\displaystyle\operatorname{BB}(\bar{\phi})=\operatorname{HCRB}(\bar{\phi})=\max_{\phi}\frac{(\phi-\bar{\phi})^{2}}{e^{\frac{2}{\sigma^{2}}\|\mathbf{m}(\phi)-\mathbf{m}(\bar{\phi})\|^{2}}-1}.$
(62)

Figure 5: 12 element uniform circular array with a radius of
$\frac{5}{3}\lambda$ and target of interest at $5\lambda$ range.

Figure 6: Near-field and far-field performance of a 12 element uniform
circular array, for a target at $5\lambda$ distance.

The results given in Figure 6 indicate that the proposed MSE expression again
predicts the threshold SNR quite closely and tracks MCRLB in the small error
region. On the other hand BB is optimistic about the threshold SNR and both
Fessler’s and So et al.’s methods yield the same results as CRLB and MCRLB for
the no model mismatch and model mismatch cases respectively.

###  7.3 DOA Estimation by an IDE (ESPRIT)

We consider a uniform linear array composed of $N=15$ sensors with $\lambda/2$
element spacing. The signal model is as follows

$\displaystyle x_{n}$ $\displaystyle=\underbrace{\alpha
e^{j\pi\cos(\bar{\phi})n}}_{\triangleq m_{n}(\bar{\phi})}+w_{n},\quad
n=0,1,\ldots,N-1,$ (63)

where $\alpha\in\mathbb{C}$ is the unknown complex amplitude,
$w_{n}\sim\mathcal{CN}(w_{n},0,\sigma_{w}^{2})$, and $\bar{\phi}=35\pi/180$
rad is the unknown true azimuth angle to be estimated. We denote the spatial
frequency with $\bar{\omega}$ and define
$\bar{\omega}\triangleq\pi\cos(\bar{\phi}).$

Due to the structure of uniform linear arrays, we can write
$m_{n}(\bar{\phi})=e^{j\bar{\omega}}m_{n-1}(\bar{\phi})$ for the elements of
the array manifold vector $m_{n}(\bar{\phi})$ in (63), which is the rotational
invariance property exploited in ESPRIT [38]. Using this property we can
define a somewhat adhoc cost function as follows

$\displaystyle J(\omega)=\sum_{n=1}^{N-1}|x_{n}-e^{j\omega}x_{n-1}|^{2}.$ (64)

By minimizing (64), we can get an estimate for $\omega$ as
$\hat{\omega}\triangleq\arg\min_{\omega}J(\omega)=\arg\big{(}\sum_{n=1}^{N-1}x_{n-1}^{*}x_{n}\big{)}$;
from which an estimate for the DOA can be generated as
$\hat{\phi}\triangleq\arccos\big{(}\frac{\hat{\omega}}{\pi}\big{)}$, which we
call the ESPRIT estimate. Note that the cost function $J(\cdot)$ is neither
symmetric around the estimate, nor is unimodal. Hence it does not satisfy the
conditions for which the proposed method yields the true MSE. The cost
function $J(\omega)$ in (64) can be written in matrix form as follows.

$J(\omega)=\left\|\mathbf{A}_{1}\mathbf{x}-e^{j\omega}\mathbf{A}_{0}\mathbf{x}\right\|^{2}=\mathbf{x}^{\mathrm{H}}{}\left(\mathbf{A}_{1}-e^{j\omega}\mathbf{A}_{0}\right)^{\mathrm{H}}{}\left(\mathbf{A}_{1}-e^{j\omega}\mathbf{A}_{0}\right)\mathbf{x},$
(65)

where

$\mathbf{x}=\left[\begin{array}[]{cccc}x_{0}&x_{1}&\ldots&x_{N-1}\end{array}\right]^{\mathrm{T}}{},\quad\mathbf{A}_{0}=\left[\begin{array}[]{cc}\mathbf{I}_{(N-1)}&\mathbf{0}_{(N-1)\times
1}\end{array}\right],\quad\mathbf{A}_{1}=\left[\begin{array}[]{cc}\mathbf{0}_{(N-1)\times
1}&\mathbf{I}_{(N-1)}\end{array}\right].$ (66)

Using $\omega\triangleq\pi\cos(\phi)$, we get,

$J(\phi)=\mathbf{x}^{\mathrm{H}}{}\big{(}\mathbf{A}_{1}-e^{j\pi\cos(\phi)}\mathbf{A}_{0}\big{)}^{\mathrm{H}}{}\big{(}\mathbf{A}_{1}-e^{j\pi\cos(\phi)}\mathbf{A}_{0}\big{)}\mathbf{x}.$
(67)

Note that in order to use the approximate MSE expression in (8), we need to
evaluate the following probability,

$\displaystyle P\left(J(\bar{\phi}+2\epsilon)\leq J(\bar{\phi})\right)$
$\displaystyle=P\left(J(\bar{\phi}+2\epsilon)-J(\bar{\phi})\leq
0\right)\triangleq P(\Delta J_{2\epsilon}\leq 0),$ (68)

where the inequalities are the reverse of those in Remark 1 since we have a
minimization problem instead of a maximization problem in our IDE. Using (67)
and after some basic algebraic operations we can express $\Delta
J_{2\epsilon}\triangleq J(\bar{\phi}+2\epsilon)-J(\bar{\phi})$ as $\Delta
J_{2\epsilon}=\mathbf{x}^{\mathrm{H}}{}\mathbf{Q}\mathbf{x}$ where

$\displaystyle\mathbf{Q}$
$\displaystyle\triangleq(e^{j\pi\cos(\bar{\phi})}-e^{j\pi\cos(\bar{\phi}+2\epsilon)})\mathbf{A}_{1}^{\mathrm{H}}{}\mathbf{A}_{0}+(e^{-j\pi\cos(\bar{\phi})}-e^{-j\pi\cos(\bar{\phi}+2\epsilon)})\mathbf{A}_{0}^{\mathrm{H}}{}\mathbf{A}_{1}.$
(69)

Even though the density of the quadratic form $\Delta
J_{2\epsilon}=\mathbf{x}^{\mathrm{H}}{}\mathbf{Q}\mathbf{x}$ is known to be
the generalized chi-squared distribution and can be evaluated numerically [39,
Appendix A], we pursue a Gaussian fit to the density in order to simplify the
probability calculations. To do that, we evaluate the first two moments of
$\Delta J_{2\epsilon}$. Using the fact that
$\mathbf{x}^{\mathrm{H}}{}\mathbf{Q}\mathbf{x}$ is always real, we can reach
the following expressions (after some algebra)

$\displaystyle\mu_{\Delta}(\bar{\phi},\epsilon)=$
$\displaystyle\,\mathbb{E}\\{\Delta
J_{2\epsilon}\\}=\,\sigma_{w}^{2}\mathrm{tr}\big{(}\mathbf{Q})+\mathbf{m}^{\mathrm{H}}{}(\bar{\phi})\mathbf{Q}\mathbf{m}(\bar{\phi}),$
(70a) $\displaystyle\sigma^{2}_{\Delta}(\bar{\phi},\epsilon)=$
$\displaystyle\,\mathrm{Var}\\{\Delta
J_{2\epsilon}\\}=\,\sigma_{w}^{4}\mathrm{tr}(\mathbf{Q}^{2})+2\sigma_{w}^{2}\mathbf{m}^{\mathrm{H}}{}(\bar{\phi})\mathbf{Q}^{2}\mathbf{m}(\bar{\phi}),$
(70b)

where
$\mathbf{m}(\phi)\triangleq[m_{0}(\phi),\,m_{1}(\phi),\cdots,m_{N-1}(\phi)]^{\mathrm{T}}{}$
and we used the result
$\mathbb{E}[(\mathbf{y}^{\mathrm{H}}{}\mathbf{Q}\mathbf{y})^{2}]=\mathrm{tr}^{2}(\mathbf{Q}\boldsymbol{\Sigma})+\mathrm{tr}((\mathbf{Q}\boldsymbol{\Sigma})^{2})$
for any Hermitian matrix $\mathbf{Q}$ and
$\mathbf{y}\sim\mathcal{CN}(\mathbf{y};\mathbf{0},\boldsymbol{\Sigma})$ [40,
Ch. V, Lemma 2.2]. With the Gaussian fit, an approximation to the suggested
MSE expression becomes

$\widehat{\operatorname{MSE}}(\bar{\phi})\approx
2\int_{-\frac{\bar{\phi}}{2}}^{\frac{\pi-\bar{\phi}}{2}}|\epsilon|\mathcal{N}_{\mathrm{cdf}}\left(0;\mu_{\Delta}(\bar{\phi},\epsilon),\sigma^{2}_{\Delta}(\bar{\phi},\epsilon)\right)\mathrm{d}\epsilon.$
(71)

where we used (9) with $\phi_{\min}=0$; $\phi_{\max}=\pi$ to set the
integration limits and the cdf of the normal distribution is used instead of
the ccdf due to the reversal of the inequalities in (68). Figure 8 shows the
results of 10,000 Monte Carlo runs for this experiment. The CRLB, BB with
single test point optimized over a grid, Fessler’s [27] and So et al.’s [28]
methods are also illustrated for comparison purposes. Note that the estimator
in this experiment is not efficient, hence its performance does not reach CRLB
at high SNR.

Figure 7: Non-random DOA estimation performance of ESPRIT along with values of
different bounds and MSE prediction expressions.

Figure 8: Bayesian DOA estimation performance of MAP and ML estimators along
with the values of BCRLB, ZZB and the proposed MSE prediction expressions (for
ML and MAP).

Consequently, the estimator performance is not characterized by the CRLB in
any SNR region. Therefore, one needs the asymptotic MSE values as well as the
pairwise error probabilities in order to calculate the MSE prediction using
MIE. BB provides a very optimistic prediction for this specific problem as in
the earlier examples. Although Fessler’s and So et al.’s methods predicted the
estimator performance well at high SNR region, they have difficulty in
representing gross errors of the estimator for low SNR values. The proposed
method, on the other hand, closely follows the estimator performance in all
SNR regions.

###  7.4 Bayesian DOA Estimation

We consider the DOA estimation problem in Section 7.3 in a Bayesian framework.
The unknown angle $\phi$ has now a prior density $f(\phi)$, which is given as
the symmetric beta distribution

$f(\phi)=\frac{1}{\pi\beta(a,a)}\left(\frac{\phi}{\pi}\right)^{a-1}\left(\frac{\pi-\phi}{\pi}\right)^{a-1},\quad\beta(a,b)\triangleq\int_{0}^{1}\phi^{a-1}(1-\phi)^{b-1}\mathrm{d}\phi,\quad
0\leq\phi\leq\pi$ (72)

with $a=10$, and the performance of the ML and MAP estimators is examined. The
proposed Bayesian MSE expression for the ML estimator can be expressed as
$\widehat{\operatorname{MSE}}_{\mathrm{ML}}=\int_{0}^{\pi}f(\phi)\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\phi)\mathrm{d}\phi$
where

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\phi)=2\int_{\frac{-\phi}{2}}^{\frac{\pi-\phi}{2}}|\epsilon|\mathcal{N}_{\mathrm{ccdf}}\left(\|\mathbf{\tilde{m}}(\phi+2\epsilon;\phi)\|;0,2\sigma_{w}^{2}\right)\mathrm{d}\epsilon,$
(73)

where the integration limits are selected as in (9) with
$\phi_{\mathrm{min}}\triangleq 0$ and $\phi_{\mathrm{max}}\triangleq\pi$. The
proposed MSE expression for the MAP estimator can be expressed as
$\widehat{\operatorname{MSE}}_{\mathrm{MAP}}=\int_{0}^{\pi}f(\phi)\widehat{\operatorname{MSE}}_{\mathrm{MAP}}(\phi)\mathrm{d}\phi$
where

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{MAP}}(\phi)=2\int_{-\pi}^{\pi}|\epsilon|\mathcal{N}_{\mathrm{ccdf}}\left(\|\mathbf{\tilde{m}}(\phi+2\epsilon;\phi)\|+\frac{\sigma_{w}^{2}}{\|\mathbf{\tilde{m}}(\phi+2\epsilon;\phi)\|}\log\left(\frac{f(\phi)}{f(\phi+2\epsilon)}\right);0,2\sigma_{w}^{2}\right)\mathrm{d}\epsilon.$
(74)

Note that $\widehat{\operatorname{MSE}}_{\mathrm{MAP}}(\phi)$ in (74) reduces
to $\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\phi)$ in (73) when the prior
is flat. BCRLB for this Bayesian estimation problem is given as [2]

$\displaystyle\mathrm{BCRLB}=\bigg{(}\pi^{2}\mathrm{SNR}\frac{N(N-1)(2N-1)}{3}\int_{0}^{\pi}\sin^{2}\phi
f(\phi)\;\mathrm{d}{}\phi+\frac{4(a-1)(2a-1)}{\pi^{2}(a-2)}\bigg{)}^{-1},$
(75)

where $\operatorname{SNR}\triangleq\frac{|\alpha|^{2}}{\sigma_{w}^{2}}$ and
$N=15$. ZZB (without the valley filling function) for the problem can be
expressed as

$\displaystyle\mathrm{ZZB}=\frac{1}{2}\int_{0}^{\pi}\int_{0}^{\pi}h(f(\phi)+f(\phi+h))P_{\mathrm{min}}^{e}(\phi,\phi+h)\mathrm{d}\theta\mathrm{d}h,$
(76)

where the minimum error probability $P_{\mathrm{min}}^{e}(\phi_{1},\phi_{2})$
can be calculated as

$\displaystyle P_{\mathrm{min}}^{e}(\phi_{1},\phi_{2})=$
$\displaystyle\,\pi_{1}\mathcal{N}_{\mathrm{ccdf}}\left(\|\mathbf{\tilde{m}}(\phi_{2};\phi_{1})\|+\frac{\sigma_{w}^{2}}{\|\mathbf{\tilde{m}}(\phi_{2};\phi_{1})\|}\log\frac{\pi_{1}}{\pi_{2}};0,2\sigma_{w}^{2}\right)$
$\displaystyle\quad+\pi_{2}\mathcal{N}_{\mathrm{ccdf}}\left(\|\mathbf{\tilde{m}}(\phi_{1};\phi_{2})\|+\frac{\sigma_{w}^{2}}{\|\mathbf{\tilde{m}}(\phi_{1};\phi_{2})\|}\log\frac{\pi_{2}}{\pi_{1}};0,2\sigma_{w}^{2}\right),$
(77)

with the prior probabilities
$\pi_{1}\triangleq\frac{f(\phi_{1})}{f(\phi_{1})+f(\phi_{2})}$ and
$\pi_{2}\triangleq 1-\pi_{1}$.

Figure 8 shows the RMSE performances of the MAP and ML estimators over 10,000
Monte Carlo runs for each $\operatorname{SNR}$ value along with the values of
BCRLB, ZZB and the proposed MSE prediction expressions
$\widehat{\operatorname{MSE}}_{\mathrm{ML}}$ and
$\widehat{\operatorname{MSE}}_{\mathrm{MAP}}$. The values of $\mathrm{ZZB}$
and the proposed $\operatorname{MSE}$ prediction expression
$\widehat{\operatorname{MSE}}_{\mathrm{MAP}}$ are identical, as expected from
the results of Section 4.3.

##  8 Conclusions

In this study we propose an MSE expression for the performance prediction of
IDEs of non-random parameters. The method provides the exact MSE value when
the objective function of the IDE is unimodal and symmetric. Even though, this
is a rather stringent restriction for the general practice; the symmetric
unimodal objective function assumption is in alignment with the operation of
consistent estimators in the asymptotic region. The maximum likelihood
estimator is the prime example for the consistent estimators. Specific to the
maximum likelihood estimator, it has been shown that the suggested MSE
expression reduces to the CRLB and MCRLB in no-misspecification and
misspecification cases, respectively. Furthermore, the suggested expression
also yields the ZZB when an a-priori distribution is assigned to the unknown
parameter for the MAP estimator.

An extension of the suggested MSE expression to the parameter estimation in
the presence of nuisance parameters is given. Numerically friendly, but
approximate, versions of the MSE expression are developed and some application
examples are given. Numerical results show that the expression not only
predicts the performance in the asymptotic region, but also provides valuable
information in the threshold region. We consider that the applicability of the
expression in other regions is related with the gradual degradation of
asymptotic region operation conditions as the operating point moves from
asymptotic region to the threshold region, say, with the reduction of SNR.

A possible interpretation for the MSE expression can be given in relation with
the method of intervals (MIE). The MIE predicts the MSE by taking into account
both small and gross error events via CRLB and the interval error
probabilities, respectively. The suggested MSE expression for the ML estimator
uses the likelihood ratio for the same purpose; but, it does not have a
problem specific interval selection.

Another interpretation for the MSE expression can be given in connection with
the ZZB. As in ZZB, the suggested MSE expression is based on the pairwise
error probabilities. Furthermore, the average of the expression for the MAP
estimator exactly reproduces ZZB for random parameters. Hence, the suggested
MSE expression for the ML estimator can also be considered, at least
informally, as the non-random parameter version of the ZZB.

An interesting observation in the non-random parameter case was that, for
medium and low SNR, the proposed MSE expressions usually slightly
underestimated the true MSE of the estimators. Hence, a potential future study
is to investigate whether the proposed expressions have any lower bounding
properties in medium and/or low SNR regions under some conditions.

##  Appendix A Proof of Theorem 1

Note that the true estimator statistic $V_{\hat{\theta}}(\theta)$ defined in
(4) can be written as [7, 41]

$\displaystyle V_{\hat{\theta}}(\theta)$
$\displaystyle=\,2\int_{0}^{\infty}\epsilon
P\big{(}|\hat{\theta}-\theta|\geq\epsilon\big{)}\;\mathrm{d}{}\epsilon=\,2\int_{0}^{\infty}\epsilon\big{[}P\big{(}\hat{\theta}-\theta\geq\epsilon\big{)}+P\big{(}\hat{\theta}-\theta\leq-\epsilon\big{)}\big{]}\;\mathrm{d}{}\epsilon,$
(78)

and the expression (5) follows from (78) if the equalities

$\displaystyle P(\hat{\theta}-\theta\geq\epsilon)=$
$\displaystyle\,P\left(\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right),$
(79a) $\displaystyle P(\hat{\theta}-\theta\leq-\epsilon)=$
$\displaystyle\,P\left(\mathcal{L}(\mathbf{x};\theta-2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right),$
(79b)

hold for $\epsilon>0$. In the following we first show that the equalities in
(79) indeed hold under the symmetry and unimodality assumptions of Theorem 1.
Only the proof of the equality $\eqref{eqn:tobeproven}$ will be made since the
proof for (79b) is very similar. In order to prove $\eqref{eqn:tobeproven}$,
we will show that $\hat{\theta}-\theta\geq\epsilon$ if and only if
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$.
The proof has two parts.

  * 1.

Proof of the implication $\hat{\theta}-\theta\geq\epsilon$ $\Rightarrow$
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$:
Suppose that $\hat{\theta}\geq\theta+\epsilon$. Since $\epsilon>0$, it is
clear that $\theta<\hat{\theta}$. If $\theta+2\epsilon<\hat{\theta}$, since
$\mathcal{L}(\mathbf{x};\theta)$ is strictly increasing for all
$\theta<\hat{\theta}$ and since $\theta<\theta+2\epsilon<\hat{\theta}$, we
would have
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)>\mathcal{L}(\mathbf{x};\theta)$ and
this would make the inequality
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$
hold. Hence, we only need to consider the case
$\theta<\hat{\theta}\leq\theta+2\epsilon$. In this case we will show that the
inequality
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$
holds by contraposition. Suppose that the reverse inequality, i.e.,
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)<\mathcal{L}(\mathbf{x};\theta)$,
holds. By the symmetry property we have

$\displaystyle\mathcal{L}(\mathbf{x};\theta+2\epsilon)=$
$\displaystyle\,\mathcal{L}(\mathbf{x};\hat{\theta}+(\theta+2\epsilon-\hat{\theta}))=\,\mathcal{L}(\mathbf{x};\hat{\theta}-(\theta+2\epsilon-\hat{\theta}))=\,\mathcal{L}(\mathbf{x};2\hat{\theta}-\theta-2\epsilon),$
(80)

which shows that
$\mathcal{L}(\mathbf{x};2\hat{\theta}-\theta-2\epsilon)<\mathcal{L}(\mathbf{x};\theta)$.
Since $2\hat{\theta}-\theta-2\epsilon\leq\hat{\theta}$ (since
$2\hat{\theta}-\theta-2\epsilon$ is the mirror image of $\theta+2\epsilon$
(with respect to $\hat{\theta}$), which is greater than or equal to
$\hat{\theta}$) and since $\mathcal{L}(\mathbf{x};\theta)$ is strictly
increasing for all $\theta<\hat{\theta}$, the inequality
$\mathcal{L}(\mathbf{x};2\hat{\theta}-\theta-2\epsilon)<\mathcal{L}(\mathbf{x};\theta)$
implies that $2\hat{\theta}-\theta-2\epsilon<\theta$. This inequality is
equivalent to the inequality $\hat{\theta}-\theta<\epsilon$, which completes
the proof.

  * 2.

Proof of the implication
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$
$\Rightarrow$ $\hat{\theta}-\theta\geq\epsilon$: Suppose that
$\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$.
Since $\mathcal{L}(\mathbf{x};\theta)$ is strictly decreasing for all
$\theta>\hat{\theta}$, we cannot have $\theta>\hat{\theta}$. Hence, we need to
have $\theta\leq\hat{\theta}$. If $\theta+2\epsilon<\hat{\theta}$, we have
$\hat{\theta}-\theta>2\epsilon>\epsilon$, which makes the inequality
$\hat{\theta}-\theta\geq\epsilon$ hold. Hence, we only need to consider the
case $\theta\leq\hat{\theta}\leq\theta+2\epsilon$. By the symmetry property
(80) we see that
$\mathcal{L}(\mathbf{x};2\hat{\theta}-\theta-2\epsilon)=\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)$.
Since we have $2\hat{\theta}-\theta-2\epsilon\leq\hat{\theta}$ and
$\theta\leq\hat{\theta}$ and since $\mathcal{L}(\mathbf{x};\theta)$ is
increasing for all $\theta<\hat{\theta}$, we need to have
$2\hat{\theta}-\theta-2\epsilon\geq\theta$. This inequality is equivalent to
the inequality $\hat{\theta}-\theta\geq\epsilon$, which completes the proof.

Hence the equalities in (79) hold and we can write (78) as

$\displaystyle V_{\hat{\theta}}(\theta)\triangleq$
$\displaystyle\,2\int_{0}^{\infty}\epsilon\big{[}P\left(\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)+P\left(\mathcal{L}(\mathbf{x};\theta-2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)\big{]}\;\mathrm{d}{}\epsilon,$
(81a) $\displaystyle=$ $\displaystyle\,2\int_{0}^{\infty}\epsilon
P\left(\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)\;\mathrm{d}{}\epsilon+2\int_{0}^{\infty}\epsilon
P\left(\mathcal{L}(\mathbf{x};\theta-2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)\;\mathrm{d}{}\epsilon,$
(81b) $\displaystyle=$ $\displaystyle\,2\int_{0}^{\infty}\epsilon
P\left(\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)\;\mathrm{d}{}\epsilon+2\int_{0}^{-\infty}\epsilon
P\left(\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)\;\mathrm{d}{}\epsilon,$
(81c) $\displaystyle=$ $\displaystyle\,2\int_{0}^{\infty}\epsilon
P\left(\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)\;\mathrm{d}{}\epsilon-2\int^{0}_{-\infty}\epsilon
P\left(\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)\;\mathrm{d}{}\epsilon,$
(81d) $\displaystyle=$
$\displaystyle\,2\int^{\infty}_{-\infty}|\epsilon|P\left(\mathcal{L}(\mathbf{x};\theta+2\epsilon)\geq\mathcal{L}(\mathbf{x};\theta)\right)\;\mathrm{d}{}\epsilon\triangleq\widehat{V}_{\hat{\theta}}(\theta),$
(81e)

which completes the proof.

##  Appendix B Proof of Proposition 1

We can write the probability in the integrand of
$\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})$ in (13) as

$\displaystyle
P\bigg{(}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
1\bigg{)}$
$\displaystyle=P\bigg{(}\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
0\bigg{)}=P\bigg{(}\frac{1}{N}\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
0\bigg{)}$ (82a)
$\displaystyle=P\bigg{(}\frac{1}{N}\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}+D(\bar{\theta}||\bar{\theta}+2\epsilon)\geq
D(\bar{\theta}||\bar{\theta}+2\epsilon)\bigg{)}$ (82b) $\displaystyle\leq
P\bigg{(}\bigg{|}\frac{1}{N}\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}+D(\bar{\theta}||\bar{\theta}+2\epsilon)\bigg{|}\geq
D(\bar{\theta}||\bar{\theta}+2\epsilon)\bigg{)}\rightarrow 0$ (82c)

for $\epsilon\neq 0$ as $N\rightarrow\infty$ where
$D(\bar{\theta}||\bar{\theta}+2\epsilon)$ stands for
$D(f(x;\bar{\theta})||f(x;\bar{\theta}+2\epsilon))$. This is because we have

$\displaystyle\frac{1}{N}\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x}|\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x}|\bar{\theta})$}}\overset{\text{p}}{\rightarrow}-D(\bar{\theta}||\bar{\theta}+2\epsilon)$
(83)

as $N\rightarrow\infty$ by the law of large numbers and
$D(\bar{\theta}||\bar{\theta}+2\epsilon)>0$ for $\epsilon\neq 0$ due to the
assumption A3. As a result, as $N\rightarrow\infty$, the integration in (13)
will be effectively only over an infinitesimal neighborhood of $\epsilon=0$
and it is only the behavior of the probability
$P\big{(}{f(\mathbf{x};\bar{\theta}+2\epsilon)}/{f(\mathbf{x};\bar{\theta})}\geq
1\big{)}$ as $\epsilon\rightarrow 0$ which determines the MSE expression
$\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})$ in (13).

Using the assumption A1, we can now obtain the Taylor expansion of $\ln
f(\mathbf{x};\bar{\theta}+2\epsilon)$ around $\epsilon=0$ given as

$\displaystyle\ln$ $\displaystyle f(\mathbf{x};\bar{\theta}+2\epsilon)=\ln
f(\mathbf{x};\bar{\theta})+2\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\bar{\theta})\epsilon+2\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\bar{\theta})\epsilon^{2}+\frac{4}{3}\frac{\partial^{3}}{\partial\theta^{3}}\ln
f(\mathbf{x};\tilde{\theta})\epsilon^{3},$ (84)

where $\tilde{\theta}$ is between $\bar{\theta}$ and $\bar{\theta}+2\epsilon$.
Since the $\frac{\partial^{3}}{\partial\theta^{3}}\ln
f(\mathbf{x};\tilde{\theta})$ is bounded by assumption A2 as
$N\rightarrow\infty$, the approximation

$\displaystyle\ln f(\mathbf{x};\bar{\theta}+2\epsilon)\approx$
$\displaystyle\ln
f(\mathbf{x};\bar{\theta})+2\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\bar{\theta})\epsilon+2\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\bar{\theta})\epsilon^{2}$ (85)

becomes valid as $\epsilon\rightarrow 0$. By rearranging, we can write

$\displaystyle\ln\frac{f(\mathbf{x};\bar{\theta}+2\epsilon)}{f(\mathbf{x};\bar{\theta})}\approx
2\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\bar{\theta})\epsilon+2\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\bar{\theta})\epsilon^{2}$ (86)

as $\epsilon\rightarrow 0$.

We can also write the Taylor expansion of $\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\hat{\theta})$ around $\theta=\bar{\theta}$ given as

$\displaystyle 0=\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\hat{\theta})=$
$\displaystyle\,\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\bar{\theta})+\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\bar{\theta})(\hat{\theta}-\bar{\theta})+\frac{1}{2}\frac{\partial^{3}}{\partial\theta^{3}}\ln
f(\mathbf{x};\theta^{\prime})(\hat{\theta}-\bar{\theta})^{2},$ (87)

where $\theta^{\prime}$ is between $\hat{\theta}$ and $\bar{\theta}$. Since
the $\frac{\partial^{3}}{\partial\theta^{3}}\ln f(\mathbf{x};\theta)$ is
bounded by assumption A2 as $N\rightarrow\infty$, the approximation

$\displaystyle 0\approx$ $\displaystyle\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\bar{\theta})+\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\bar{\theta})(\hat{\theta}-\bar{\theta})$ (88)

becomes valid as $\hat{\theta}\overset{\text{a.s.}}{\rightarrow}\bar{\theta}$
as $N\rightarrow\infty$. Rearranging, we obtain

$\displaystyle\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\bar{\theta})\approx-\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\bar{\theta})(\hat{\theta}-\bar{\theta})$ (89)

as $N\rightarrow\infty$. Substituting $\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\bar{\theta})$ in (89) into (86), we get

$\displaystyle\ln\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\approx-2\epsilon\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\bar{\theta})(\hat{\theta}-\bar{\theta}-\epsilon)$ (90)

as $\epsilon\rightarrow 0$ and $N\rightarrow\infty$. We can now substitute the
result (90) into the the probability in the integrand of
$\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})$ in (13) to obtain

$\displaystyle P$
$\displaystyle\bigg{(}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
1\bigg{)}=P\bigg{(}-\epsilon\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\bar{\theta})(\hat{\theta}-\bar{\theta}-\epsilon)\geq 0\bigg{)}.$
(91)

Using assumptions A1-A4, it can be shown that (See [33, Lemma 2.1 Part-i] or
[34, Lemma 4.1 Part-i])

$\displaystyle-\frac{1}{N}\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\bar{\theta})\overset{\text{p}}{\rightarrow}$
$\displaystyle-\mathbb{E}\bigg{[}\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(x;\bar{\theta})\bigg{]}=\,\mathcal{I}(\bar{\theta})>0$ (92)

as $N\rightarrow\infty$ where we used the law of large numbers. This allows us
to write (91) as

$\displaystyle
P\bigg{(}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
1\bigg{)}=$
$\displaystyle\,P\big{(}\epsilon(\hat{\theta}-\bar{\theta}-\epsilon)\geq
0\big{)}=\begin{cases}P\big{(}\hat{\theta}-\bar{\theta}\geq\epsilon\big{)},&\epsilon>0\\\
1,&\epsilon=0\\\
P\big{(}\hat{\theta}-\bar{\theta}\leq\epsilon\big{)},&\epsilon<0\end{cases},$
(93)

for $\epsilon\rightarrow 0$ and $N\rightarrow\infty$. Note that the
probabilities $P\big{(}\hat{\theta}-\bar{\theta}\geq\epsilon\big{)}$,
$\epsilon>0$ and $P\big{(}\hat{\theta}-\bar{\theta}\leq\epsilon\big{)}$,
$\epsilon<0$ would vanish as $N\rightarrow\infty$, just as the probability
$P\big{(}f(\mathbf{x};\bar{\theta}+2\epsilon)/f(\mathbf{x};\bar{\theta})\geq
1\big{)}$, $\epsilon\neq 0$, itself, thanks to the fact that
$\hat{\theta}\overset{\text{a.s.}}{\rightarrow}\bar{\theta}$. As a result, we
can substitute the right hand side of (93) into the finite support version of
the integral (13) to get

$\displaystyle\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\bar{\theta})$
$\displaystyle\triangleq
2\int_{\frac{\theta_{\mathrm{min}}-\bar{\theta}}{2}}^{\frac{\theta_{\mathrm{max}}-\bar{\theta}}{2}}|\epsilon|P\bigg{(}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
1\bigg{)}\;\mathrm{d}{}\epsilon=2\int_{\theta_{\mathrm{min}}-\bar{\theta}}^{\theta_{\mathrm{max}}-\bar{\theta}}|\epsilon|P\bigg{(}\genfrac{}{}{}{}{\raisebox{1.5pt}{$f(\mathbf{x};\bar{\theta}+2\epsilon)$}}{\raisebox{-1.5pt}{$f(\mathbf{x};\bar{\theta})$}}\geq
1\bigg{)}\;\mathrm{d}{}\epsilon$ (94a)
$\displaystyle=-2\int_{\theta_{\mathrm{min}}-\bar{\theta}}^{0}\epsilon
P\big{(}\hat{\theta}-\bar{\theta}\leq\epsilon\big{)}\;\mathrm{d}{}\epsilon+2\int_{0}^{\theta_{\mathrm{max}}-\bar{\theta}}\epsilon
P\big{(}\hat{\theta}-\bar{\theta}\geq\epsilon\big{)}\;\mathrm{d}{}\epsilon$
(94b)
$\displaystyle=\int_{\theta_{\mathrm{min}}-\bar{\theta}}^{0}\epsilon^{2}f_{\hat{\theta}-\bar{\theta}}(\epsilon)\;\mathrm{d}{}\epsilon+\int_{0}^{\theta_{\mathrm{max}}-\bar{\theta}}\epsilon^{2}f_{\hat{\theta}-\bar{\theta}}(\epsilon)\;\mathrm{d}{}\epsilon=\int_{\theta_{\mathrm{min}}-\bar{\theta}}^{\theta_{\mathrm{max}}-\bar{\theta}}\epsilon^{2}f_{\hat{\theta}-\bar{\theta}}(\epsilon)\;\mathrm{d}{}\epsilon=\mathbb{E}[(\hat{\theta}-\bar{\theta})^{2}]\triangleq\operatorname{MSE}_{\mathrm{ML}}(\bar{\theta}),$
(94c)

as $N\rightarrow\infty$, which completes the proof of (17). The proof of (18)
follows trivially if the ML estimate $\hat{\theta}$ is also asymptotically
efficient.

##  Appendix C Proof of Proposition 2

In the case of MML estimation, we set
$\mathcal{L}(\mathbf{x};\theta)\triangleq\ln f(\mathbf{x};\theta)$ in the
probability in the integrand of $\widehat{V}_{\hat{\theta}}(\theta)$ in (6).
We can now write the probability in $\widehat{V}_{\hat{\theta}}(\theta)$ as

$\displaystyle P\big{(}$ $\displaystyle\ln
f(\mathbf{x};\theta_{*}+2\epsilon)\geq\ln
f(\mathbf{x};\theta_{*})\big{)}=P\bigg{(}\frac{1}{N}\ln\frac{f(\mathbf{x};\theta_{*}+2\epsilon)}{f(\mathbf{x};\theta_{*})}\geq
0\bigg{)}$ (95a) $\displaystyle=$
$\displaystyle\,P\bigg{(}\frac{1}{N}\ln\frac{\bar{f}(\mathbf{x})}{f(\mathbf{x};\theta_{*})}-\frac{1}{N}\ln\frac{\bar{f}(\mathbf{x})}{f(\mathbf{x};\theta_{*}+2\epsilon)}\geq
0\bigg{)}$ (95b) $\displaystyle=$
$\displaystyle\,P\bigg{(}\frac{1}{N}\ln\frac{\bar{f}(\mathbf{x})}{f(\mathbf{x};\theta_{*})}-\frac{1}{N}\ln\frac{\bar{f}(\mathbf{x})}{f(\mathbf{x};\theta_{*}+2\epsilon)}-\big{(}D(\theta_{*})-D(\theta_{*}+2\epsilon)\big{)}\geq
D(\theta_{*}+2\epsilon)-D(\theta_{*})\bigg{)}$ (95c) $\displaystyle\leq$
$\displaystyle\,P\bigg{(}\bigg{|}\frac{1}{N}\ln\frac{\bar{f}(\mathbf{x})}{f(\mathbf{x};\theta_{*})}-\frac{1}{N}\ln\frac{\bar{f}(\mathbf{x})}{f(\mathbf{x};\theta_{*}+2\epsilon)}-\big{(}D(\theta_{*})-D(\theta_{*}+2\epsilon)\big{)}\bigg{|}\geq
D(\theta_{*}+2\epsilon)-D(\theta_{*})\bigg{)}\rightarrow 0$ (95d)

for $\epsilon\neq 0$ as $N\rightarrow\infty$ where $D(\theta)$ stands for
$D(\bar{f}(x)||f(x;\theta))$. This is because we have

$\displaystyle\frac{1}{N}\ln\frac{\bar{f}(\mathbf{x})}{f(\mathbf{x};\theta_{*})}-\frac{1}{N}\ln\frac{\bar{f}(\mathbf{x})}{f(\mathbf{x};\theta_{*}+2\epsilon)}\overset{\text{p}}{\rightarrow}D(\theta_{*})-D(\theta_{*}+2\epsilon)$
(96)

as $N\rightarrow\infty$ by the law of large numbers and
$D(\theta_{*}+2\epsilon)>D(\theta_{*})$ for $\epsilon\neq 0$ due to the
assumption A3. As a result, as $N\rightarrow\infty$, the integration in (8)
will be effectively only over an infinitesimal neighborhood of $\epsilon=0$
and it is only the behavior of the probability $P\big{(}\ln
f(\mathbf{x};\theta_{*}+2\epsilon)\geq\ln f(\mathbf{x};\theta_{*})\big{)}$ as
$\epsilon\rightarrow 0$ which determines the expression
$\widehat{V}_{\hat{\theta}}(\theta_{*})$.

Following a similar approach that is used for obtaining (86), we can write

$\displaystyle\ln\frac{f(\mathbf{x};\theta_{*}+2\epsilon)}{f(\mathbf{x};\theta_{*})}\approx
2\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\theta_{*})\epsilon+2\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\theta_{*})\epsilon^{2}$ (97)

as $\epsilon\rightarrow 0$. Using an approach similar to that used for
obtaining (89) we can get

$\displaystyle\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\theta_{*})\approx-\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\theta_{*})(\hat{\theta}-\theta_{*})$ (98)

as $N\rightarrow\infty$. Substituting $\frac{\partial}{\partial\theta}\ln
f(\mathbf{x};\theta_{*})$ in (98) into (97), we get

$\displaystyle\ln\frac{f(\mathbf{x};\theta_{*}+2\epsilon)}{f(\mathbf{x};\theta_{*})}\approx-2\epsilon\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\theta_{*})(\hat{\theta}-\theta_{*}-\epsilon)$ (99)

as $\epsilon\rightarrow 0$ and $N\rightarrow\infty$. We can now substitute the
result (99) into the probability in the integrand of $\widehat{V}(\theta_{*})$

$\displaystyle P$ $\displaystyle\big{(}\ln
f(\mathbf{x};\theta_{*}+2\epsilon)\geq\ln
f(\mathbf{x};\theta_{*})\big{)}=P\bigg{(}-\epsilon\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\theta_{*})(\hat{\theta}-\theta_{*}-\epsilon)\geq 0\bigg{)}.$
(100)

Using the assumptions A1-A4, it can be shown that (See [33, Lemma 2.1 Part-i]
or [34, Lemma 4.1 Part-i])

$\displaystyle-\frac{1}{N}\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(\mathbf{x};\theta_{*})\overset{\text{p}}{\rightarrow}$
$\displaystyle-\mathbb{E}_{\bar{f}}\bigg{[}\frac{\partial^{2}}{\partial\theta^{2}}\ln
f(x;\theta_{*})\bigg{]}=-\mathcal{A}(\theta_{*})>0$ (101)

as $N\rightarrow\infty$ where we used the law of large numbers. This allows us
to write (100) as

$\displaystyle P\big{(}\ln f(\mathbf{x};\theta_{*}+2\epsilon)\geq$
$\displaystyle\,\ln
f(\mathbf{x};\theta_{*})\big{)}=\,P\big{(}\epsilon(\hat{\theta}-\theta_{*}-\epsilon)\geq
0\big{)}=\begin{cases}P\big{(}\hat{\theta}-\theta_{*}\geq\epsilon\big{)},&\epsilon>0\\\
1,&\epsilon=0\\\
P\big{(}\hat{\theta}-\theta_{*}\leq\epsilon\big{)},&\epsilon<0\end{cases},$
(102)

for $\epsilon\rightarrow 0$ and $N\rightarrow\infty$. Note that the
probabilities $P\big{(}\hat{\theta}-\theta_{*}\geq\epsilon\big{)}$,
$\epsilon>0$ and $P\big{(}\hat{\theta}-\theta_{*}\leq\epsilon\big{)}$,
$\epsilon<0$ would vanish as $N\rightarrow\infty$, just as the probability
$P\big{(}\ln f(\mathbf{x};\theta_{*}+2\epsilon)\geq\ln
f(\mathbf{x};\theta_{*})\big{)}$, $\epsilon\neq 0$, itself, thanks to the fact
that $\hat{\theta}\overset{\text{a.s.}}{\rightarrow}\theta_{*}$. As a result,
we can substitute the right hand side of (102) into the finite support version
of the integral in $\widehat{V}_{\hat{\theta}}(\theta_{*})$ to get

$\displaystyle\widehat{V}_{\text{MML}}(\theta_{*})\triangleq$ $\displaystyle
2\int_{\frac{\theta_{\mathrm{min}}-\theta_{*}}{2}}^{\frac{\theta_{\mathrm{max}}-\theta_{*}}{2}}|\epsilon|P\left(\ln
f(\mathbf{x};\theta_{*}+2\epsilon)\geq\ln
f(\mathbf{x};\theta_{*})\right)\;\mathrm{d}{}\epsilon$ (103a)
$\displaystyle=2\int_{\theta_{\mathrm{min}}-\theta_{*}}^{\theta_{\mathrm{max}}-\theta_{*}}|\epsilon|P\left(\ln
f(\mathbf{x};\theta_{*}+2\epsilon)\geq\ln
f(\mathbf{x};\theta_{*})\right)\;\mathrm{d}{}\epsilon$ (103b)
$\displaystyle=-2\int_{\theta_{\mathrm{min}}-\theta_{*}}^{0}\epsilon
P\big{(}\hat{\theta}-\theta_{*}\leq\epsilon\big{)}\;\mathrm{d}{}\epsilon+2\int_{0}^{\theta_{\mathrm{max}}-\theta_{*}}\epsilon
P\big{(}\hat{\theta}-\theta_{*}\geq\epsilon\big{)}\;\mathrm{d}{}\epsilon$
(103c)
$\displaystyle=\int_{\theta_{\mathrm{min}}-\theta_{*}}^{0}\epsilon^{2}f_{\hat{\theta}-\theta_{*}}(\epsilon)\;\mathrm{d}{}\epsilon+\int_{0}^{\theta_{\mathrm{max}}-\theta_{*}}\epsilon^{2}f_{\hat{\theta}-\theta_{*}}(\epsilon)\;\mathrm{d}{}\epsilon$
(103d)
$\displaystyle=\int_{\theta_{\mathrm{min}}-\theta_{*}}^{\theta_{\mathrm{max}}-\theta_{*}}\epsilon^{2}f_{\hat{\theta}-\bar{\theta}}(\epsilon)\;\mathrm{d}{}\epsilon$
(103e) $\displaystyle=\mathbb{E}[(\hat{\theta}-\theta_{*})^{2}]\triangleq
V_{\text{MML}}(\theta_{*}),$ (103f)

as $N\rightarrow\infty$, which completes the proof of (20). The proof of (21)
follows trivially if the MML estimate $\hat{\theta}$ is also asymptotically
misspecified efficient.

##  Appendix D Implementation Details of the Methods Used in Section 7

###  D.1 Implementation Details for Section 7.1

  * 1.

ML estimate is calculated on a grid. The grid consists of 3600 uniformly
spaced points in the interval $[-\pi,\,\pi]$ for (scalar) azimuth estimation,
and 3600 uniformly spaced points in the interval $[0,\,\pi]$ for (scalar)
elevation estimation. For the case of nuisance parameter, azimuth and
elevation angles are selected such that the corresponding unit vectors have an
approximately uniform distribution over the unit sphere. To do so, elevation
interval $[0,\,\pi]$ is divided into 200 equally spaced values
$\\{\theta_{k}\\}_{k=1}^{200}$, and each constant elevation circle on the unit
sphere is divided into $\lceil{(100\sin(\theta_{k})}\rceil$ points where
$\lceil{\cdot}\rceil$ denotes the ceiling function. By doing so, a total of
$198192$ grid points ($\phi$-$\theta$ pairs) are obtained.

BBs [14] are calculated with a single test point which is optimized on a grid.
The single and multiple parameter BB grids are selected the same as the grids
used for ML. MIE grid points are selected as the local maxima outside the
mainlobe in the beampattern for the scalar parameter estimation case, which
requires a peak finding algorithm over the 1D beampattern function. For the
case of a nuisance parameter, the beampattern function is a 2D surface, and
grid points are the local maxima on this surface. For both cases the
beampattern function is calculated over the ML grid, and local maxima are
found using Matlab built-in functions (findpeaks($\cdot$) for maxima on the 1D
curve, imregionalmax($\cdot$) for maxima on the 2D surface).

Proposed method does not require grid point selection for the scalar parameter
case. For the case of a nuisance parameter, grid points are selected over the
entire support of the nuisance parameter with logarithmic spacing around the
true value, making the grid denser as the grid points approach the true
parameter value. In Matlab notation, the grid is defined as follows

$\displaystyle\boldsymbol{\delta\theta}$
$\displaystyle=\texttt{logspace}(-7,\mathrm{log_{10}}e_{\text{max}},60),$
(104a) $\displaystyle\boldsymbol{\Theta}_{\backslash 1}$
$\displaystyle=[\bar{\theta},\,\bar{\theta}-\boldsymbol{\delta\theta},\,\bar{\theta}+\boldsymbol{\delta\theta}].$
(104b)

The statement (104a) generates 60 logarithmically spaced points within
$[10^{-7},e_{\text{max}}]$ interval in Matlab. We selected
$e_{\text{max}}=\pi/2$ for the elevation angle as the nuisance parameter and
$e_{\text{max}}=\pi$ for the azimuth angle as the nuisance parameter.

###  D.2 Implementation Details for Section 7.2

  * 1.

ML grid consists of 3600 uniformly spaced points in the interval
$[-\pi,\,\pi]$ for (scalar) azimuth estimation. BB grid is the same as the ML
grid. Other methods do not require grid points.

###  D.3 Implementation Details for Section 7.3

  * 1.

BB grid is the same as the ML grid for the problem in Section 7.2. Other
methods do not require grid points.

The numerical integrals of the proposed method in all subsections up to this
point are taken using the Matlab function integral($\cdot$) with the following
tolerance values: AbsTol=1e-5, RelTol=1e-5. The integration limits for the
numerical integrals of the proposed method were set as
$\epsilon\in\big{[}\frac{-\pi-\bar{\phi}}{2},\frac{\pi-\bar{\phi}}{2}\big{]}$
for azimuth estimation and
$\epsilon\in\big{[}-\frac{\bar{\theta}}{2},\frac{\pi-\bar{\theta}}{2}\big{]}$
for elevation estimation.

###  D.4 Implementation Details for Section 7.4

  * 1.

ML and MAP use the same grid, which consists of 8192 uniformly spaced points
over $[-\pi,\pi]$ for $\omega=\pi\cos(\phi)$. The (non-uniform) grid points
for $\phi$ are calculated using expression $\arccos(\frac{\omega}{\pi})$ from
the grid for $\omega$.

  * 2.

The double integral for ZZB is taken using the MATLAB function
integral2($\cdot$) with the default tolerance settings.

  * 3.

The numerical integrals of $\widehat{\operatorname{MSE}}_{\mathrm{ML}}(\phi)$
and $\widehat{\operatorname{MSE}}_{\mathrm{MAP}}(\phi)$ are taken using the
Matlab function integral($\cdot$) with the tolerance values AbsTol=1e-18,
RelTol=1e-12. The integrals with respect to the prior are calculated over a
uniform grid over the inverval $[0,\pi]$ with the grid spacing 0.01.

## References

  * [1] S. M. Kay, Fundamentals of Statistical Signal Processing, Vol I: Estimation Theory, Prentice-Hall, 1993.
  * [2] H. L. V. Trees, K. L. Bell, Bayesian Bounds for Parameter Estimation and Nonlinear Filtering/Tracking, Wiley-IEEE Press, 2007.
  * [3] H. L. V. Trees, Detection, Estimation and Modulation Theory, John Wiley and Sons Inc., 1968.
  * [4] A. J. Weiss, E. Weinstein, Fundamental limitations in passive time delay estimation - Part I : Narrow-band systems, IEEE Trans. Acoust. Speech Signal Process. 31 (2) (1983) 472–486.
  * [5] J. Ziv, M. Zakai, Some lower bounds on signal parameter estimation, IEEE Trans. Inf. Theory 15 (3) (1969) 386–391.
  * [6] S. Bellini, G. Tartara, Bounds on error in signal parameter estimation, IEEE Trans. Comm. 22 (1974) 340–342.
  * [7] K. L. Bell, Y. Steinberg, Y. Ephraim, H. L. Van Trees, Extended Ziv-Zakai lower bound for vector parameter estimation, IEEE Trans. Inf. Theory 43 (2) (1997) 624–637.
  * [8] E. Nitzan, T. Routtenberg, J. Tabrikian, A new class of Bayesian cyclic bounds for periodic parameter estimation, IEEE Trans. Sig. Proc. 64 (1) (2016) 229–243.
  * [9] A. Xu, T. Coleman, Minimax lower bounds for circular source localization, in: 2020 IEEE Int. Symp. on Inf. Theory (ISIT), 2020, pp. 1242–1247.
  * [10] T. Routtenberg, J. Tabrikian, Bayesian periodic Cramér-Rao bound, IEEE Sig. Proc. Lett. 29 (2022) 1878–1882.
  * [11] R. G. Gallager, Stochastic Processes: Theory for Applications, Cambridge Uni. Press, New York, 2014.
  * [12] D. G. Chapman, H. Robbins, Minimum variance estimation without regularity assumptions, Ann. Math. Statist. 22 (4) (1951) 581–586.
  * [13] E. W. Barankin, Locally best unbiased estimates, Ann. Math. Statist. 20 (1949) 477–501.
  * [14] R. McAulay, E. Hofstetter, Barankin bounds on parameter estimation, IEEE Trans. Inf. Theory 17 (6) (1971) 669–676.
  * [15] P. Forster, P. Larzabal, On lower bounds for deterministic parameter estimation, in: IEEE Int. Conf. on Acoustics, Speech, and Sig. Proc., Vol. 2, 2002, pp. II–1137–II–1140.
  * [16] K. Todros, J. Tabrikian, General classes of performance lower bounds for parameter estimation - Part I: Non-Bayesian bounds for unbiased estimators, IEEE Trans. Inf. Theory 56 (10) (October 2010).
  * [17] T. Routtenberg, J. Tabrikian, Non-Bayesian periodic Cramér-Rao bound, IEEE Trans. Sig. Proc. 61 (4) (2013) 1019–1032.
  * [18] T. Routtenberg, J. Tabrikian, Cyclic Barankin-type bounds for non-Bayesian periodic parameter estimation, IEEE Trans. Sig. Proc. 62 (13) (2014) 3321–3336.
  * [19] T. Routtenberg, J. Tabrikian, Cyclic Cramér-Rao-type bounds for periodic parameter estimation, in: 2016 19th Int. Conf. on Inf. Fusion (FUSION), 2016, pp. 1797–1804.
  * [20] K. Todros, R. Winik, J. Tabrikian, On the limitations of Barankin type bounds for MLE threshold prediction, Signal Processing 108 (2015) 622–627.
  * [21] P. Stoica, R. L. Moses, On biased estimators and the unbiased Cramer-Rao lower bound, Signal Processing 21 (4) (1990) 349–350.
  * [22] H. White, Consequences and detection of misspecified nonlinear regression models, Journal of the American Statistical Association 76 (374) (1981) 419–433.
  * [23] L. Wasserman, All of Statistics: A Concise Course in Statistical Inference, Springer, New York, 2010.
  * [24] S. Fortunati, F. Gini, M. S. Greco, C. D. Richmond, Performance bounds for parameter estimation under misspecified models: Fundamental findings and applications, IEEE Sig. Proc. Mag. 34 (6) (2017) 142–157.
  * [25] C. D. Richmond, L. L. Horowitz, Parameter bounds on estimation accuracy under model misspecification, IEEE Trans. Sig. Proc. 63 (9) (2015) 2263–2278.
  * [26] H. White, Maximum likelihood estimation of misspecified models, Econometrica 50 (1) (1982) 1–25.
  * [27] J. A. Fessler, Mean and variance of implicitly defined biased estimators (such as penalized maximum likelihood): Applications to tomography, IEEE Trans. Image Proc. 5 (3) (1996) 493–506.
  * [28] H. C. So, Y. T. Chan, K. C. Ho, Y. Chan, Simple formulas for bias and mean square error computation, IEEE Sig. Proc. Mag. (2013) 162–165.
  * [29] F. Athley, Threshold region performance of maximum likelihood direction of arrival estimators, IEEE Trans. Sig. Proc. 53 (4) (2005) 1359–1373.
  * [30] C. D. Richmond, Mean-squared error and threshold SNR prediction of maximum likelihood signal parameter estimation with estimated colored noise covariances, IEEE Trans. Inf. Theory 52 (5) (2006) 2146–2164.
  * [31] A. Mallat, S. Gezici, D. Dardari, C. Craeye, L. Vandendorpe, Statistics of the MLE and approximate upper and lower bounds - Part I: Application to TOA estimation, IEEE Trans. Sig. Proc. 62 (21) (2014) 5663–5676.
  * [32] A. Mallat, S. Gezici, D. Dardari, C. Craeye, L. Vandendorpe, Statistics of the MLE and approximate upper and lower bounds - Part II: Threshold computation and optimal pulse design for TOA estimation, IEEE Trans. Sig. Proc. 62 (21) (2014) 5677–5689.
  * [33] Q. H. Vuong, Cramer-Rao bounds for misspecified models, https://resolver.caltech.edu/CaltechAUTHORS:20170823-162930200, unpublished (Oct. 1986).
  * [34] S. Fortunati, F. Gini, M. S. Greco, Parameter bounds under misspecified models for adaptive radar detection, Vol. 7 of Academic Press Library in Signal Processing, Academic Press, 2017, Ch. 4, pp. 197–252.
  * [35] S. M. Kay, Fundamentals of Statistical Signal Processing, Volume II: Detection Theory, Prentice-Hall, 1998.
  * [36] D. Vu, A. Renaux, R. Boyer, S. Marcos, A Cramer Rao bounds based analysis of 3D antenna array geometries made from ULA branches, Multidimensional Systems and Sig. Proc. 24 (2013) 121–155.
  * [37] H. L. V. Trees, Optimum Array Processing, Wiley and Sons Inc., 2002.
  * [38] R. Roy, T. Kailath, ESPRIT–Estimation of signal parameters via rotational invariance techniques, IEEE Trans. Acoust., Speech, Signal Proc. 37 (7) (1989) 984–995.
  * [39] M. S. Paolella, Linear Models and Time‐Series Analysis: Regression, ANOVA, ARMA and GARCH, John Wiley & Sons, 2018.
  * [40] K. S. Miller, Complex Stochastic Processes: An Introduction to Theory and Application, Addison-Wesley Publishing Company, 1974.
  * [41] K. L. Bell, Performance bounds in parameter estimation with application to bearing estimation, Ph.D. thesis, George Mason University (1995).