# Bayesian Estimation Applied to Stochastic Localization with Constraints due to Interfaces and Boundaries

Authors: Wolfgang Hoegele; Rainer Loeschel; Barbara Dobler; Oliver Koelbl; Piotr Zygmanski
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Purpose. We present a systematic Bayesian formulation of the stochastic localization/triangulation problem close to constraining interfaces. Methods. For this purpose, the terminology of Bayesian estimation is summarized suitably for applied researchers including the presentation of Maximum Likelihood (ML), Maximum A Posteriori (MAP), and Minimum Mean Square Error (MMSE) estimation. Explicit estimators for triangulation are presented for the linear 2D parallel beam and the nonlinear 3D cone beam model. The priors in MAP and MMSE optionally incorporate (A) the hard constraints about the interface and (B) knowledge about the probability of the object with respect to the interface. All presented estimators are compared in several simulation studies for live acquisition scenarios with 10,000 samples each.
Results. First, the presented application shows that MAP and MMSE perform considerably better, leading to lower Root Mean Square Errors (RMSEs) in the simulation studies compared to the ML approach by typically introducing a bias. Second, utilizing priors including (A) and (B) is very beneficial compared to just including (A). Third, typically MMSE leads to better results than MAP, by the cost of significantly higher computational effort. Conclusion. Depending on the specific application and prior knowledge, MAP and MMSE estimators strongly increase the estimation accuracy for localization close to interfaces.


---

## Body

## 1. Introduction

Due to their inherent ability to use prior knowledge, Bayesian approaches have gained interest in many fields in the recent years, such as in computer vision [1, 2] and medical applications [3–6], as well as robotics (Kalman-Filter, Particle-Filter) [7, 8] and metrology [9]. In this work, we present a Bayesian framework for the localization of objects with constraints derived from interfaces (surfaces) or boundaries specific to a given localization problem.Methods to find the locations of objects based on several projections are known from photogrammetry for a long time. These methods are also called triangulation in a multiple view geometry (e.g., see Hartley and Zisserman [1] for an overview of established methods). Introducing probabilities in order to deal with the Gaussian measurement noise is a well-known approach for estimating a 3-dimensional (3D) point based on two 2-dimensional (2D) projections, for example, see Hartley and Sturm [10]. In this context, Bedekar  and Haralick [11] provided a Maximum A Posteriori (MAP) estimator by introducing prior knowledge independent of the observation for the estimation based on a set of 2D projections, which represents a starting point for this Paper but does not provide a more general perspective on estimation which is required to bring it to the level described in this paper. We are following this probabilistic interpretation of the inverse problem of triangulation and present explicit solutions in order to characterize the triangulation problem as well as embed it into the general framework of Bayesian estimation. In particular, we demonstrate the benefit of including soft and hard constraints about the solution space into the prior of the Bayesian triangulation problem. We compare three stochastic estimation methods: Maximum Likelihood (ML), Maximum A Posteriori (MAP), and Minimum Mean Square Error (MMSE), and show that, depending on the prior knowledge, MAP or MMSE is a suitable approach for this kind of problem. We show how the estimation methods can be utilized for a 3D nonlinear/non-Gaussian geometrical model in physical space without modifications or approximations of the observation geometry.In the localization problem the input data consists of 2D projections (or images) of the 3D objects and the goal is the optimal localization of features contained in these objects, by employing prior knowledge. In general, projections of 3D objects onto 2D planes contain information about their coordinates, but due to the experimental conditions these coordinates may not be accurate or may not be easily extractable from these projections. Thus the localization problem is typically an ill-posed inverse problem. Fortunately, in many applications there exists prior knowledge about the probability densities of these uncertainties. Furthermore, if there is a general probabilistic knowledge about the position of the observed object or additional constraints about the domain of its coordinates, the total prior knowledge becomes difficult to be incorporated into the localization model. Such localizations close to a constraining interface or boundary can be modeled and solved in the Bayesian approach. For instance, in the field of medicine, during image-guided procedures (interventional radiology, radiotherapy, and surgical procedures with C-arm CBXT scanners) there is often a need to obtain 3D coordinates of anatomical features or implanted markers instantly from multiple X-ray projections, which are constrained by external or internal anatomical interfaces (e.g., patient surface) and form a hard or semi-flexible barrier. Both, the general Bayesian framework and its application to a specific localization problem presented herein can also be utilized for other applications in image or data processing, such as computer vision or metrology.The general aim of this paper is to provide a pragmatic foundation for the Bayesian estimation in a compact way suitable for applied sciences utilizing a concise and uniform mathematical notation (general parts about the Bayesian methodology of this paper are partly taken from the doctoral thesis of the first author, which is submitted (and not yet published) at the Regensburg University, Germany, at the time of this paper submission), which is difficult to find in more theoretical formulations and whose full understanding is required for the localization problem at stake. This brief introduction to Bayesian estimation may benefit several research fields in which there is a gap between the general mathematical concepts and direct applications of Bayesian methods. Thus the major motivation of presenting the Bayesian estimation theory is the linking of these two areas by providing the Bayesian formalism for researches in applied sciences. The second aim of this work is to demonstrate that a Bayesian framework can be effectively utilized for more accurate localization of features close to interfaces, which may have many applications.The structure of the paper is as follows: first, we present the basic terminology of Bayesian estimation required to derive the solutions presented in the subsequent sections based on a concise summary of the books of Lehmann and Casella [12] and Robert [13]. Second, we specifically present the application of stochastic localization applied to 1D and 2D projections of 2D and 3D objects. Besides such a structured presentation of Bayesian triangulation, the hope is that this course of explanation will help the readers to develop their own Bayesian estimators for their specific problems at hand.


## 2. General Bayesian Theory

### 2.1. Introduction to the Terminology of Estimation

The goal of estimation theory is to determine a concrete set of unknown real-valued parametersx∈ΩX⊆ℝm (e.g., the 2D/3D position coordinates of the observed object in triangulation) from real-valued observations u∈ΩU⊆ℝn (e.g., the detected or measured 1D/2D projections of the object), which are imperfectly or ambiguously related to each other (bold fonts of symbols are utilized to emphasize multidimensional vectors). Their ambiguous relation is due to the experimental uncertainties of the given problem, such as an imperfect Radon transform, or similar transforms. The estimation of parameters x from observables u demands a set of definitions which will be briefly introduced in this section. The concrete parameters x and observations u are associated with continuous  random variablesX and U, respectively (denoted by corresponding capital letters), which are further associated with the underlying probability distributions describing the sources of uncertainties. With this concept a connection between the deterministic and probabilistic quantities is established. For instance, the probability that x takes a value lower than s is defined by
(1)P(X≤s)≔∫-∞s‍f[X](x)dx,
with the corresponding Probability Density Function (PDF) f[X](x). This fact is summarized by stating that X~f[X].In general, the occurrence of one event can influence the probability of another. It is possible to define such a residual probability by  conditional probabilities, and for the corresponding PDF describing the random variableU under the condition of X=x this is formulated by
(2)f[U](u∣x)≔f[U,X](u,x)f[X](x).
In this formula, the expression f[U,X](u,x) describes the joint PDF of U and X. Based on this concept, the most important quantities for estimation are defined in the next subsections.


#### 2.1.1. Likelihood, Prior, and Posterior

In  Table1 the most important types of PDFs for estimation are presented and interpreted in the following paragraphs. The PDF f[U](u∣x) defines the probability of the current observations u if the true values x are known. In a typical observation scenario this corresponds to the accuracy of the observation device with respect to a specific true parameter x. It is common to define the  likelihood function ℒu(x)≔f[U](u∣x) describing the same probabilities of the observation but this time with the focus on the conditioning parameter x. Therefore the likelihood is regarded as a similarity measure between the current observation u and possibly “true” parameters x.Table 1
Probability Density Functions (PDF) which are important for estimation theory.
f
[
U
]
(
u
∣
x
)
PDF of the observationu conditioned on x
ℒ
u
(
x
)
:
=
f
[
U
]
(
u
∣
x
)
The likelihood function (depending on the condition x)
π
(
x
)
:
=
f
[
X
]
(
x
)
theprior, that is, PDF of the prior knowledge about x
π
(
x
∣
u
)
:
=
f
[
X
]
(
x
∣
u
)
theposterior, that is, PDF of the posterior distribution of parameter x conditioned on u
∫
Ω
X
f
[
U
]
(
u
∣
x
)
·
π
(
x
)
d
x
PDF of the probability of the occurrence ofu (=: marginal distribution)
X
~
f
[
X
]
,
U
~
f
[
U
]
random variablesX and U are distributed according to their PDFs f[X] and f[U]The PDF of thepriorπ(x)≔f[X](x) defines how probable the occurrence of x is; that is, X~π. The use of the prior is an intensely discussed topic, since one could derive arbitrary conclusions by utilizing a specific and strong prior knowledge which supports a certain point of view which may be just one of many possible options. Although this is generally valid, the utilization of appropriate prior knowledge (noninformative priors, priors based on actual observations of the same or similar situations, or priors derived from knowledge about the range of x) can reduce the uncertainties inherent to x in addition to making observations u.InBayesian statistics the most important PDF is the  posteriorπ(x∣u)≔f[X](x∣u), which defines the probabilities of x under the condition of the observation u. Similar to the likelihood this quantity is also often regarded as a similarity measure between the observation and x.The aforementioned PDFs are related to each other byBayes’ Theorem, which is expressed as
(3)π(x∣u)=f[U](u∣x)·π(x)∫ΩX‍f[U](u∣x)·π(x)dx.
The benefit of this theorem is that it allows to reverse the condition by introducing prior knowledge about the parameter of interest. In practical use, with this theorem the posteriorπ(x∣u) can be determined by having only knowledge about the likelihoodℒu(x) and the priorπ(x), since the denominator on the right side is typically determined indirectly by normalizing the posterior to the integral value of 1.


#### 2.1.2. Loss and Risk

Anestimator is any function δ:ΩU→ΩX, which utilizes the current observations u in order to find an approximation to the true parameter x. In order to decide on the use of an estimator compared to another, one needs a quality criterion, which depends on the true value for the parameters x and an estimate δ(u) based on the observation u. In estimation theory this is realized with the concept of the loss functionL(x,δ(u)). An important subset is strictly convex loss functions, such as the popular quadratic loss,
(4)L(x,δ(u))≔∥x-δ(u)∥2,
since they lead to unique Bayesian estimators. An example for a not strictly but still convex loss is the constant loss(5)L(x,δ(u))≔{0,formaxi∈{1,…,p}|xi-δi(u)|≤ε,1,else,
for a constant ε>0, which is of technical interest for the derivation of the MAP estimator. In this context, a very attractive thought is to find an optimal estimator δ(·) by minimizing the loss function L(x,δ(u)). Unfortunately, the naive minimization of the loss would lead to δ(u)=x, which is a practically useless estimator since it demands knowledge about the true x. Essentially, the loss itself does not contain any information about the relation between the parameter x and the measurements u and, therefore, represents an incomplete description of the estimation problem. Due to this observation, in classical estimation theory the concept of the risk function is introduced, which is defined here by the conditional expectation value.Definition 1.
The  risk function  (or alternatively average loss) is defined as (It is inherent to the classical frequentist’s approach that the parameters are not interpreted as random variablesX. Since in this work a general Bayesian perspective is focused on, in the definition of the frequentist risk, the parameters are written as random variables but simultaneously conditioned on X=x, which essentially leads to the same risk function.)(6)R(x,δ)≔𝔼[L(X,δ(U))∣X=x]=∫ΩU‍L(x,δ(u))·f[U](u∣x)du.The risk functionR(x,δ) incorporates all possible outcomes u of the observation and weights them with their likelihood f[U](u∣x). For this reason, the long-term average loss of a given estimator δ(·) is obtained by the risk function of (6). In consequence, based on the risk function a reasonable comparison of estimators can be performed using the concept of admissibility (see Lehmann and Casella [12, pp. 309–376] or Robert [13, pp. 74–77, 391–415] for details). On the other side, the risk function still depends on the unknown parameter x, and, therefore, direct minimization of the risk function can only be performed for a fixed x. In order to incorporate all possible values of x into a single quality criterion an extension of the concept of the risk function is introduced by incorporating (possibly noninformative) prior knowledge about x by π(x).Definition 2.
The  integrated risk function is defined as the expectation value taken overX and U(7)r(π,δ)≔𝔼[L(X,δ(U))]=∫ΩX‍∫ΩU‍L(x,δ(u))·f[U](u∣x)·π(x)dudx.With this criterion one gets a single real value for a given estimatorδ(·) which incorporates all parameters x and possible observations u and allows a straight comparison between different estimators (and especially a meaningful optimization). Unfortunately, to find an optimal estimator directly from the minimization of the integrated risk of  (7) is a nontrivial problem but can be solved by Bayesian estimation as presented in Section 2.2. Further, a delicate point is the introduction of a prior distribution π(x) before any measurement is performed since this could bias the estimation result in a disadvantageous way. On the other side, if  ΩX is bounded and π(x)=1/vol(ΩX)=const. for all x∈ΩX, all possible values x are treated equally likely in order to erase the subjective notion when introducing the prior. In consequence, the application of the integrated risk function in order to derive estimators could also be of interest in cases in which one wants to avoid such a subjective notion and still get an estimator working best for all possible values x and u in average.


### 2.2. Maximum Likelihood, Maximum A Posteriori, and Bayesian Estimation

In this section, three central estimation methods are introduced in a general way (for a brief characterization of these approaches, see  Section2.3).


#### 2.2.1. Maximum Likelihood (ML)

Maximum Likelihood estimation is typically applied in situations in which no prior knowledge aboutx∈ΩX is available.Definition 3.
A Maximum Likelihood (ML) estimator for a given observationu∈ΩU is
(8)δ(u)≔argmaxxℒu(x)=argmaxxf[U](u∣x),
where “arg max” means that we want to find the position x which maximizes the likelihood ℒu(x). This is also valid for the logarithm of the likelihood since the logarithm is a strictly increasing function, and the likelihood is positive. Since the ML estimation depends on the shape of the likelihood function which is problem-specific, under fairly general assumptions, the existence but not uniqueness of the estimate can be guaranteed. The theorem by Lehmann and Casella [12, p. 463] assures the existence of an asymptotically unbiased, consistent, and asymptotically efficient estimator which is the solution of ∇xℒu(x)=0, the so-called likelihood equations. If there is a unique solution to these equations, then in consequence, the corresponding ML fulfills these properties. From a practical perspective the determination of the estimate results in an optimization problem, and numerically this may lead to problems with respect to running into local maxima instead of global maxima.


#### 2.2.2. Bayesian Estimator (BE)

Bayesian estimation utilizes prior knowledge aboutx∈ΩX, and focuses on the given observation. For this reason, a modified risk function needs to be introduced at first.Definition 4.
The  Bayesian risk function is defined for a given lossL, observation u∈ΩU and the posterior π(x∣u) as
(9)ρu(π,δ)≔𝔼[L(X,δ(U))∣U=u]=∫ΩX‍L(x,δ(u))·π(x∣u)dx.The main difference between the Bayesian risk (9) and the integrated risk (7) is that the integrated risk does average over all possible observations u∈ΩU, and the Bayesian risk utilizes only the concretely acquired observation u.Definition 5.
A Bayesian estimator (BE) for a given observationu∈ΩU is the estimator δ(u) which minimizes the Bayesian risk function
(10)δ(u)≔argminξρu(π,ξ).At the first glance, the Bayesian estimator therefore looks unrelated to the optimization of the previously introduced risk functions in  Section2.1. A powerful statement is achieved by the following theorem.Theorem 6.
The estimation function which minimizes the integrated riskr(π,δ) can be obtained by selecting for every specific observation u∈ΩU the estimate δ(u) which minimizes the Bayesian risk ρu(π,δ).For the proof see Robert [13, p. 63]. A major property of Bayesian estimators is that if the loss is strictly convex, then the corresponding BE is unique. This will be shown explicitly in the next paragraph for the quadratic loss. For this loss the properties are similar as for the ML estimator, leading to consistent and asymptotically efficient estimators under some weak assumptions on the prior; see Lehmann and Casella [12, p. 490]. From a practical point of view, the numerical computation of the integral might be limiting the applications of BE if analytical solutions cannot be derived. Compared to the optimization problem in ML, the general calculation of these multiple integrals is typically computationally more challenging.


#### 2.2.3. The Minimum Mean Square Error (MMSE)

The Bayesian estimator utilizing the quadratic loss is called the Minimum Mean Square Error (MMSE). The derivation leads to(11)δ(u)≔argminξ𝔼[L(X,ξ)∣U=u]=argminξ∫ΩX‍∥x-ξ∥2·π(x∣u)dx.
Since the loss is strictly convex and  (11) represents a parameter integral with respect to ξ, the unique solution can be derived from
(12)∇ξ𝔼[L(X,ξ)∣U=u]=!0.
This means for every coordinate i∈{1,…,m}(13)∫ΩX‍2(xi-ξi)·π(x∣u)dx=!0,ξi·∫ΩX‍π(x∣u)dx︸=1=∫ΩX‍xi·π(x∣u)dx.
This leads to the corresponding MMSE for each coordinate i∈{1,…,m} with
(14)δi(u)=𝔼[Xi∣U=u]=c·∫ΩX‍xi·ℒu(x)·π(x)dx,
which is the expectation value of the posterior for every coordinate i∈{1,…,m}, with the normalization constant c. This example already points out the importance of the shape of the posterior π(x∣u) for the Bayesian approach.


#### 2.2.4. Maximum A Posteriori (MAP)

An alternative approach, conceptually in between ML and BE, is the Maximum A Posteriori (MAP) estimation which already utilizes prior knowledge aboutx∈ΩX, but has many characteristics of the ML approach.Definition 7.
The Maximum A Posteriori estimator (MAP) for a given observationu∈ΩU is
(15)δ(u)≔argmaxxπ(x∣u)=argmaxxℒu(x)·π(x),
since the marginal distribution is positive and independent of x. In other words, MAP is the mode of the posterior.When applying the logarithm toℒu(x)·π(x), this approach can be identified with a specific Tikhonov regularization of the ML approach (also known as penalized Maximum Likelihood), as discussed in Eldar [14]. Since the MAP estimation depends on the shape of the likelihood and prior function which is problem specifically defined, under fairly general assumptions, the existence but not uniqueness of the estimate can be guaranteed. The asymptotic properties are the same as for the ML case, with some further regularity conditions on ℒu(x) and π(x); for example, see Robert [13, p. 166]. Although MAP utilizes prior knowledge about x∈ΩX, the numerical problem is the same as for ML, that is, resulting in a numerical optimization problem with the problem of identifying local and global maxima.In the context of Bayesian estimation, the MAP estimation can be derived as a limit of Bayesian estimators with respect to the constant loss function  (5). It is shown that
(16)δ(u)≔argminξ𝔼[L(X,ξ)∣U=u]=argminξ∫ΩX‍L(x,ξ)·π(x∣u)dx=argmaxξ∫Cξ‍π(x∣u)dx,
with Cξ as the m-dimensional cube with side length 2ε and midpoint ξ. In other words, the aim is to find the midpoint of the m-dimensional cube Cξ, so that the volume under the graph of the posterior is maximized. By assuming an at least continuous posterior and letting ε→0, the best midpoint δ(u) is the mode of the posterior, that is, the MAP estimator. This is illustrated in  Figure 1. Since the constant loss function is not strictly convex, the corresponding BEs are not necessarily unique, and this implicates the nonuniqueness of the MAP as a limit of these BEs.Figure 1
The MAP estimator as a limit of Bayesian estimators withx∈ℝ2. For every fixed ε the Bayesian estimator is the midpoint δ(u) which maximizes the gray shaded volume.


### 2.3. Characterization of the Estimators

The presented estimators are closely related which is graphically illustrated in  Figure2.  A major distinction is drawn in between the use of prior knowledge about x∈ΩX or not, which is shown as a thick boundary in  Figure 2. In detail, the main connections are summarized in the following list.  (i)
BE→ MMSE: the general BE reduces to the MMSE when applying the quadratic loss of  (4).(ii)
MMSE→ MAP: if the posterior distribution is symmetric, that is, π(μ+x∣u)=π(μ-x∣u) around the global maximum μ, then the MMSE estimator is identical with the MAP estimator. In this case, due to symmetry the expectation and the maximum of the posterior coincide.(iii)
BE→ MAP: the BE converges to the MAP estimator when applying the constant loss of  (5) with ε→0.(iv)
MAP→ ML: if the prior is noninformative, that is, constant for a finite ΩX or asymptotically noninformative for an infinite ΩX as the deviation parameters of the prior go to infinity, then MAP coincides with the ML estimate.Figure 2
Illustration of the relation between estimation based on ML, MAP, MMSE, and Bayesian estimation. The intersections state under which conditions the estimators are equivalent. Every condition appears only once; that is, since the conditions are not excluding each other, multiple conditions lead to further equalities.Of course, the combination of the stated properties is possible, for example, for a symmetric posterior distribution around a global maximum and a noninformative prior the ML and MMSE coincide. There are several other concepts utilized in the literature for deriving non-Bayesian estimators, that is, without (possibly noninformative) prior knowledge aboutx∈ΩX, such as the following. (i)
Unbiasedness, that is, the search of an unbiased estimator which minimizes the frequentist risk locally or uniformly. This property has been regarded as very important from the beginning of estimation theory (e.g., see Cramér [15]), although it could be shown that an unbiased estimator does not always exist (see, e.g., Lehmann and Casella [12, p. 83]), and lower risk function values (6) can be achieved when introducing a suitable bias (Eldar [16]). In this context, the Minimum Variance Unbiased Estimator (MVUE) is of special interest, since it represents the efficient unbiased estimator.(ii)
Minimaxity, that is, the minimization of the maximum value of the risk function of  (6), which typically leads to rather conservative estimates, since they can be approximated by Bayesian estimators with a prior distributions concentrated at the worst case (Robert [13, p. 67]). Unfortunately, these estimates do not always exist (Robert [13, p. 69]).(iii)
Equivariance, that is, the invariance of the estimation problem under certain transformations of the estimator. This property is important for cases in which a given estimator is applied to modified estimation situations, but the subsequent investigation of the optimality of the modified estimator should be avoided. This is important for practical cases, in which the estimator is modified routinely or a certain symmetry of the estimator is demanded by the nature of the parameters one wants to estimate. For details, see Lehmann and Casella [12, p. 147–225].


## 2.1. Introduction to the Terminology of Estimation

The goal of estimation theory is to determine a concrete set of unknown real-valued parametersx∈ΩX⊆ℝm (e.g., the 2D/3D position coordinates of the observed object in triangulation) from real-valued observations u∈ΩU⊆ℝn (e.g., the detected or measured 1D/2D projections of the object), which are imperfectly or ambiguously related to each other (bold fonts of symbols are utilized to emphasize multidimensional vectors). Their ambiguous relation is due to the experimental uncertainties of the given problem, such as an imperfect Radon transform, or similar transforms. The estimation of parameters x from observables u demands a set of definitions which will be briefly introduced in this section. The concrete parameters x and observations u are associated with continuous  random variablesX and U, respectively (denoted by corresponding capital letters), which are further associated with the underlying probability distributions describing the sources of uncertainties. With this concept a connection between the deterministic and probabilistic quantities is established. For instance, the probability that x takes a value lower than s is defined by
(1)P(X≤s)≔∫-∞s‍f[X](x)dx,
with the corresponding Probability Density Function (PDF) f[X](x). This fact is summarized by stating that X~f[X].In general, the occurrence of one event can influence the probability of another. It is possible to define such a residual probability by  conditional probabilities, and for the corresponding PDF describing the random variableU under the condition of X=x this is formulated by
(2)f[U](u∣x)≔f[U,X](u,x)f[X](x).
In this formula, the expression f[U,X](u,x) describes the joint PDF of U and X. Based on this concept, the most important quantities for estimation are defined in the next subsections.


### 2.1.1. Likelihood, Prior, and Posterior

In  Table1 the most important types of PDFs for estimation are presented and interpreted in the following paragraphs. The PDF f[U](u∣x) defines the probability of the current observations u if the true values x are known. In a typical observation scenario this corresponds to the accuracy of the observation device with respect to a specific true parameter x. It is common to define the  likelihood function ℒu(x)≔f[U](u∣x) describing the same probabilities of the observation but this time with the focus on the conditioning parameter x. Therefore the likelihood is regarded as a similarity measure between the current observation u and possibly “true” parameters x.Table 1
Probability Density Functions (PDF) which are important for estimation theory.
f
[
U
]
(
u
∣
x
)
PDF of the observationu conditioned on x
ℒ
u
(
x
)
:
=
f
[
U
]
(
u
∣
x
)
The likelihood function (depending on the condition x)
π
(
x
)
:
=
f
[
X
]
(
x
)
theprior, that is, PDF of the prior knowledge about x
π
(
x
∣
u
)
:
=
f
[
X
]
(
x
∣
u
)
theposterior, that is, PDF of the posterior distribution of parameter x conditioned on u
∫
Ω
X
f
[
U
]
(
u
∣
x
)
·
π
(
x
)
d
x
PDF of the probability of the occurrence ofu (=: marginal distribution)
X
~
f
[
X
]
,
U
~
f
[
U
]
random variablesX and U are distributed according to their PDFs f[X] and f[U]The PDF of thepriorπ(x)≔f[X](x) defines how probable the occurrence of x is; that is, X~π. The use of the prior is an intensely discussed topic, since one could derive arbitrary conclusions by utilizing a specific and strong prior knowledge which supports a certain point of view which may be just one of many possible options. Although this is generally valid, the utilization of appropriate prior knowledge (noninformative priors, priors based on actual observations of the same or similar situations, or priors derived from knowledge about the range of x) can reduce the uncertainties inherent to x in addition to making observations u.InBayesian statistics the most important PDF is the  posteriorπ(x∣u)≔f[X](x∣u), which defines the probabilities of x under the condition of the observation u. Similar to the likelihood this quantity is also often regarded as a similarity measure between the observation and x.The aforementioned PDFs are related to each other byBayes’ Theorem, which is expressed as
(3)π(x∣u)=f[U](u∣x)·π(x)∫ΩX‍f[U](u∣x)·π(x)dx.
The benefit of this theorem is that it allows to reverse the condition by introducing prior knowledge about the parameter of interest. In practical use, with this theorem the posteriorπ(x∣u) can be determined by having only knowledge about the likelihoodℒu(x) and the priorπ(x), since the denominator on the right side is typically determined indirectly by normalizing the posterior to the integral value of 1.


### 2.1.2. Loss and Risk

Anestimator is any function δ:ΩU→ΩX, which utilizes the current observations u in order to find an approximation to the true parameter x. In order to decide on the use of an estimator compared to another, one needs a quality criterion, which depends on the true value for the parameters x and an estimate δ(u) based on the observation u. In estimation theory this is realized with the concept of the loss functionL(x,δ(u)). An important subset is strictly convex loss functions, such as the popular quadratic loss,
(4)L(x,δ(u))≔∥x-δ(u)∥2,
since they lead to unique Bayesian estimators. An example for a not strictly but still convex loss is the constant loss(5)L(x,δ(u))≔{0,formaxi∈{1,…,p}|xi-δi(u)|≤ε,1,else,
for a constant ε>0, which is of technical interest for the derivation of the MAP estimator. In this context, a very attractive thought is to find an optimal estimator δ(·) by minimizing the loss function L(x,δ(u)). Unfortunately, the naive minimization of the loss would lead to δ(u)=x, which is a practically useless estimator since it demands knowledge about the true x. Essentially, the loss itself does not contain any information about the relation between the parameter x and the measurements u and, therefore, represents an incomplete description of the estimation problem. Due to this observation, in classical estimation theory the concept of the risk function is introduced, which is defined here by the conditional expectation value.Definition 1.
The  risk function  (or alternatively average loss) is defined as (It is inherent to the classical frequentist’s approach that the parameters are not interpreted as random variablesX. Since in this work a general Bayesian perspective is focused on, in the definition of the frequentist risk, the parameters are written as random variables but simultaneously conditioned on X=x, which essentially leads to the same risk function.)(6)R(x,δ)≔𝔼[L(X,δ(U))∣X=x]=∫ΩU‍L(x,δ(u))·f[U](u∣x)du.The risk functionR(x,δ) incorporates all possible outcomes u of the observation and weights them with their likelihood f[U](u∣x). For this reason, the long-term average loss of a given estimator δ(·) is obtained by the risk function of (6). In consequence, based on the risk function a reasonable comparison of estimators can be performed using the concept of admissibility (see Lehmann and Casella [12, pp. 309–376] or Robert [13, pp. 74–77, 391–415] for details). On the other side, the risk function still depends on the unknown parameter x, and, therefore, direct minimization of the risk function can only be performed for a fixed x. In order to incorporate all possible values of x into a single quality criterion an extension of the concept of the risk function is introduced by incorporating (possibly noninformative) prior knowledge about x by π(x).Definition 2.
The  integrated risk function is defined as the expectation value taken overX and U(7)r(π,δ)≔𝔼[L(X,δ(U))]=∫ΩX‍∫ΩU‍L(x,δ(u))·f[U](u∣x)·π(x)dudx.With this criterion one gets a single real value for a given estimatorδ(·) which incorporates all parameters x and possible observations u and allows a straight comparison between different estimators (and especially a meaningful optimization). Unfortunately, to find an optimal estimator directly from the minimization of the integrated risk of  (7) is a nontrivial problem but can be solved by Bayesian estimation as presented in Section 2.2. Further, a delicate point is the introduction of a prior distribution π(x) before any measurement is performed since this could bias the estimation result in a disadvantageous way. On the other side, if  ΩX is bounded and π(x)=1/vol(ΩX)=const. for all x∈ΩX, all possible values x are treated equally likely in order to erase the subjective notion when introducing the prior. In consequence, the application of the integrated risk function in order to derive estimators could also be of interest in cases in which one wants to avoid such a subjective notion and still get an estimator working best for all possible values x and u in average.


## 2.1.1. Likelihood, Prior, and Posterior

In  Table1 the most important types of PDFs for estimation are presented and interpreted in the following paragraphs. The PDF f[U](u∣x) defines the probability of the current observations u if the true values x are known. In a typical observation scenario this corresponds to the accuracy of the observation device with respect to a specific true parameter x. It is common to define the  likelihood function ℒu(x)≔f[U](u∣x) describing the same probabilities of the observation but this time with the focus on the conditioning parameter x. Therefore the likelihood is regarded as a similarity measure between the current observation u and possibly “true” parameters x.Table 1
Probability Density Functions (PDF) which are important for estimation theory.
f
[
U
]
(
u
∣
x
)
PDF of the observationu conditioned on x
ℒ
u
(
x
)
:
=
f
[
U
]
(
u
∣
x
)
The likelihood function (depending on the condition x)
π
(
x
)
:
=
f
[
X
]
(
x
)
theprior, that is, PDF of the prior knowledge about x
π
(
x
∣
u
)
:
=
f
[
X
]
(
x
∣
u
)
theposterior, that is, PDF of the posterior distribution of parameter x conditioned on u
∫
Ω
X
f
[
U
]
(
u
∣
x
)
·
π
(
x
)
d
x
PDF of the probability of the occurrence ofu (=: marginal distribution)
X
~
f
[
X
]
,
U
~
f
[
U
]
random variablesX and U are distributed according to their PDFs f[X] and f[U]The PDF of thepriorπ(x)≔f[X](x) defines how probable the occurrence of x is; that is, X~π. The use of the prior is an intensely discussed topic, since one could derive arbitrary conclusions by utilizing a specific and strong prior knowledge which supports a certain point of view which may be just one of many possible options. Although this is generally valid, the utilization of appropriate prior knowledge (noninformative priors, priors based on actual observations of the same or similar situations, or priors derived from knowledge about the range of x) can reduce the uncertainties inherent to x in addition to making observations u.InBayesian statistics the most important PDF is the  posteriorπ(x∣u)≔f[X](x∣u), which defines the probabilities of x under the condition of the observation u. Similar to the likelihood this quantity is also often regarded as a similarity measure between the observation and x.The aforementioned PDFs are related to each other byBayes’ Theorem, which is expressed as
(3)π(x∣u)=f[U](u∣x)·π(x)∫ΩX‍f[U](u∣x)·π(x)dx.
The benefit of this theorem is that it allows to reverse the condition by introducing prior knowledge about the parameter of interest. In practical use, with this theorem the posteriorπ(x∣u) can be determined by having only knowledge about the likelihoodℒu(x) and the priorπ(x), since the denominator on the right side is typically determined indirectly by normalizing the posterior to the integral value of 1.


## 2.1.2. Loss and Risk

Anestimator is any function δ:ΩU→ΩX, which utilizes the current observations u in order to find an approximation to the true parameter x. In order to decide on the use of an estimator compared to another, one needs a quality criterion, which depends on the true value for the parameters x and an estimate δ(u) based on the observation u. In estimation theory this is realized with the concept of the loss functionL(x,δ(u)). An important subset is strictly convex loss functions, such as the popular quadratic loss,
(4)L(x,δ(u))≔∥x-δ(u)∥2,
since they lead to unique Bayesian estimators. An example for a not strictly but still convex loss is the constant loss(5)L(x,δ(u))≔{0,formaxi∈{1,…,p}|xi-δi(u)|≤ε,1,else,
for a constant ε>0, which is of technical interest for the derivation of the MAP estimator. In this context, a very attractive thought is to find an optimal estimator δ(·) by minimizing the loss function L(x,δ(u)). Unfortunately, the naive minimization of the loss would lead to δ(u)=x, which is a practically useless estimator since it demands knowledge about the true x. Essentially, the loss itself does not contain any information about the relation between the parameter x and the measurements u and, therefore, represents an incomplete description of the estimation problem. Due to this observation, in classical estimation theory the concept of the risk function is introduced, which is defined here by the conditional expectation value.Definition 1.
The  risk function  (or alternatively average loss) is defined as (It is inherent to the classical frequentist’s approach that the parameters are not interpreted as random variablesX. Since in this work a general Bayesian perspective is focused on, in the definition of the frequentist risk, the parameters are written as random variables but simultaneously conditioned on X=x, which essentially leads to the same risk function.)(6)R(x,δ)≔𝔼[L(X,δ(U))∣X=x]=∫ΩU‍L(x,δ(u))·f[U](u∣x)du.The risk functionR(x,δ) incorporates all possible outcomes u of the observation and weights them with their likelihood f[U](u∣x). For this reason, the long-term average loss of a given estimator δ(·) is obtained by the risk function of (6). In consequence, based on the risk function a reasonable comparison of estimators can be performed using the concept of admissibility (see Lehmann and Casella [12, pp. 309–376] or Robert [13, pp. 74–77, 391–415] for details). On the other side, the risk function still depends on the unknown parameter x, and, therefore, direct minimization of the risk function can only be performed for a fixed x. In order to incorporate all possible values of x into a single quality criterion an extension of the concept of the risk function is introduced by incorporating (possibly noninformative) prior knowledge about x by π(x).Definition 2.
The  integrated risk function is defined as the expectation value taken overX and U(7)r(π,δ)≔𝔼[L(X,δ(U))]=∫ΩX‍∫ΩU‍L(x,δ(u))·f[U](u∣x)·π(x)dudx.With this criterion one gets a single real value for a given estimatorδ(·) which incorporates all parameters x and possible observations u and allows a straight comparison between different estimators (and especially a meaningful optimization). Unfortunately, to find an optimal estimator directly from the minimization of the integrated risk of  (7) is a nontrivial problem but can be solved by Bayesian estimation as presented in Section 2.2. Further, a delicate point is the introduction of a prior distribution π(x) before any measurement is performed since this could bias the estimation result in a disadvantageous way. On the other side, if  ΩX is bounded and π(x)=1/vol(ΩX)=const. for all x∈ΩX, all possible values x are treated equally likely in order to erase the subjective notion when introducing the prior. In consequence, the application of the integrated risk function in order to derive estimators could also be of interest in cases in which one wants to avoid such a subjective notion and still get an estimator working best for all possible values x and u in average.


## 2.2. Maximum Likelihood, Maximum A Posteriori, and Bayesian Estimation

In this section, three central estimation methods are introduced in a general way (for a brief characterization of these approaches, see  Section2.3).


### 2.2.1. Maximum Likelihood (ML)

Maximum Likelihood estimation is typically applied in situations in which no prior knowledge aboutx∈ΩX is available.Definition 3.
A Maximum Likelihood (ML) estimator for a given observationu∈ΩU is
(8)δ(u)≔argmaxxℒu(x)=argmaxxf[U](u∣x),
where “arg max” means that we want to find the position x which maximizes the likelihood ℒu(x). This is also valid for the logarithm of the likelihood since the logarithm is a strictly increasing function, and the likelihood is positive. Since the ML estimation depends on the shape of the likelihood function which is problem-specific, under fairly general assumptions, the existence but not uniqueness of the estimate can be guaranteed. The theorem by Lehmann and Casella [12, p. 463] assures the existence of an asymptotically unbiased, consistent, and asymptotically efficient estimator which is the solution of ∇xℒu(x)=0, the so-called likelihood equations. If there is a unique solution to these equations, then in consequence, the corresponding ML fulfills these properties. From a practical perspective the determination of the estimate results in an optimization problem, and numerically this may lead to problems with respect to running into local maxima instead of global maxima.


### 2.2.2. Bayesian Estimator (BE)

Bayesian estimation utilizes prior knowledge aboutx∈ΩX, and focuses on the given observation. For this reason, a modified risk function needs to be introduced at first.Definition 4.
The  Bayesian risk function is defined for a given lossL, observation u∈ΩU and the posterior π(x∣u) as
(9)ρu(π,δ)≔𝔼[L(X,δ(U))∣U=u]=∫ΩX‍L(x,δ(u))·π(x∣u)dx.The main difference between the Bayesian risk (9) and the integrated risk (7) is that the integrated risk does average over all possible observations u∈ΩU, and the Bayesian risk utilizes only the concretely acquired observation u.Definition 5.
A Bayesian estimator (BE) for a given observationu∈ΩU is the estimator δ(u) which minimizes the Bayesian risk function
(10)δ(u)≔argminξρu(π,ξ).At the first glance, the Bayesian estimator therefore looks unrelated to the optimization of the previously introduced risk functions in  Section2.1. A powerful statement is achieved by the following theorem.Theorem 6.
The estimation function which minimizes the integrated riskr(π,δ) can be obtained by selecting for every specific observation u∈ΩU the estimate δ(u) which minimizes the Bayesian risk ρu(π,δ).For the proof see Robert [13, p. 63]. A major property of Bayesian estimators is that if the loss is strictly convex, then the corresponding BE is unique. This will be shown explicitly in the next paragraph for the quadratic loss. For this loss the properties are similar as for the ML estimator, leading to consistent and asymptotically efficient estimators under some weak assumptions on the prior; see Lehmann and Casella [12, p. 490]. From a practical point of view, the numerical computation of the integral might be limiting the applications of BE if analytical solutions cannot be derived. Compared to the optimization problem in ML, the general calculation of these multiple integrals is typically computationally more challenging.


### 2.2.3. The Minimum Mean Square Error (MMSE)

The Bayesian estimator utilizing the quadratic loss is called the Minimum Mean Square Error (MMSE). The derivation leads to(11)δ(u)≔argminξ𝔼[L(X,ξ)∣U=u]=argminξ∫ΩX‍∥x-ξ∥2·π(x∣u)dx.
Since the loss is strictly convex and  (11) represents a parameter integral with respect to ξ, the unique solution can be derived from
(12)∇ξ𝔼[L(X,ξ)∣U=u]=!0.
This means for every coordinate i∈{1,…,m}(13)∫ΩX‍2(xi-ξi)·π(x∣u)dx=!0,ξi·∫ΩX‍π(x∣u)dx︸=1=∫ΩX‍xi·π(x∣u)dx.
This leads to the corresponding MMSE for each coordinate i∈{1,…,m} with
(14)δi(u)=𝔼[Xi∣U=u]=c·∫ΩX‍xi·ℒu(x)·π(x)dx,
which is the expectation value of the posterior for every coordinate i∈{1,…,m}, with the normalization constant c. This example already points out the importance of the shape of the posterior π(x∣u) for the Bayesian approach.


### 2.2.4. Maximum A Posteriori (MAP)

An alternative approach, conceptually in between ML and BE, is the Maximum A Posteriori (MAP) estimation which already utilizes prior knowledge aboutx∈ΩX, but has many characteristics of the ML approach.Definition 7.
The Maximum A Posteriori estimator (MAP) for a given observationu∈ΩU is
(15)δ(u)≔argmaxxπ(x∣u)=argmaxxℒu(x)·π(x),
since the marginal distribution is positive and independent of x. In other words, MAP is the mode of the posterior.When applying the logarithm toℒu(x)·π(x), this approach can be identified with a specific Tikhonov regularization of the ML approach (also known as penalized Maximum Likelihood), as discussed in Eldar [14]. Since the MAP estimation depends on the shape of the likelihood and prior function which is problem specifically defined, under fairly general assumptions, the existence but not uniqueness of the estimate can be guaranteed. The asymptotic properties are the same as for the ML case, with some further regularity conditions on ℒu(x) and π(x); for example, see Robert [13, p. 166]. Although MAP utilizes prior knowledge about x∈ΩX, the numerical problem is the same as for ML, that is, resulting in a numerical optimization problem with the problem of identifying local and global maxima.In the context of Bayesian estimation, the MAP estimation can be derived as a limit of Bayesian estimators with respect to the constant loss function  (5). It is shown that
(16)δ(u)≔argminξ𝔼[L(X,ξ)∣U=u]=argminξ∫ΩX‍L(x,ξ)·π(x∣u)dx=argmaxξ∫Cξ‍π(x∣u)dx,
with Cξ as the m-dimensional cube with side length 2ε and midpoint ξ. In other words, the aim is to find the midpoint of the m-dimensional cube Cξ, so that the volume under the graph of the posterior is maximized. By assuming an at least continuous posterior and letting ε→0, the best midpoint δ(u) is the mode of the posterior, that is, the MAP estimator. This is illustrated in  Figure 1. Since the constant loss function is not strictly convex, the corresponding BEs are not necessarily unique, and this implicates the nonuniqueness of the MAP as a limit of these BEs.Figure 1
The MAP estimator as a limit of Bayesian estimators withx∈ℝ2. For every fixed ε the Bayesian estimator is the midpoint δ(u) which maximizes the gray shaded volume.


## 2.2.1. Maximum Likelihood (ML)

Maximum Likelihood estimation is typically applied in situations in which no prior knowledge aboutx∈ΩX is available.Definition 3.
A Maximum Likelihood (ML) estimator for a given observationu∈ΩU is
(8)δ(u)≔argmaxxℒu(x)=argmaxxf[U](u∣x),
where “arg max” means that we want to find the position x which maximizes the likelihood ℒu(x). This is also valid for the logarithm of the likelihood since the logarithm is a strictly increasing function, and the likelihood is positive. Since the ML estimation depends on the shape of the likelihood function which is problem-specific, under fairly general assumptions, the existence but not uniqueness of the estimate can be guaranteed. The theorem by Lehmann and Casella [12, p. 463] assures the existence of an asymptotically unbiased, consistent, and asymptotically efficient estimator which is the solution of ∇xℒu(x)=0, the so-called likelihood equations. If there is a unique solution to these equations, then in consequence, the corresponding ML fulfills these properties. From a practical perspective the determination of the estimate results in an optimization problem, and numerically this may lead to problems with respect to running into local maxima instead of global maxima.


## 2.2.2. Bayesian Estimator (BE)

Bayesian estimation utilizes prior knowledge aboutx∈ΩX, and focuses on the given observation. For this reason, a modified risk function needs to be introduced at first.Definition 4.
The  Bayesian risk function is defined for a given lossL, observation u∈ΩU and the posterior π(x∣u) as
(9)ρu(π,δ)≔𝔼[L(X,δ(U))∣U=u]=∫ΩX‍L(x,δ(u))·π(x∣u)dx.The main difference between the Bayesian risk (9) and the integrated risk (7) is that the integrated risk does average over all possible observations u∈ΩU, and the Bayesian risk utilizes only the concretely acquired observation u.Definition 5.
A Bayesian estimator (BE) for a given observationu∈ΩU is the estimator δ(u) which minimizes the Bayesian risk function
(10)δ(u)≔argminξρu(π,ξ).At the first glance, the Bayesian estimator therefore looks unrelated to the optimization of the previously introduced risk functions in  Section2.1. A powerful statement is achieved by the following theorem.Theorem 6.
The estimation function which minimizes the integrated riskr(π,δ) can be obtained by selecting for every specific observation u∈ΩU the estimate δ(u) which minimizes the Bayesian risk ρu(π,δ).For the proof see Robert [13, p. 63]. A major property of Bayesian estimators is that if the loss is strictly convex, then the corresponding BE is unique. This will be shown explicitly in the next paragraph for the quadratic loss. For this loss the properties are similar as for the ML estimator, leading to consistent and asymptotically efficient estimators under some weak assumptions on the prior; see Lehmann and Casella [12, p. 490]. From a practical point of view, the numerical computation of the integral might be limiting the applications of BE if analytical solutions cannot be derived. Compared to the optimization problem in ML, the general calculation of these multiple integrals is typically computationally more challenging.


## 2.2.3. The Minimum Mean Square Error (MMSE)

The Bayesian estimator utilizing the quadratic loss is called the Minimum Mean Square Error (MMSE). The derivation leads to(11)δ(u)≔argminξ𝔼[L(X,ξ)∣U=u]=argminξ∫ΩX‍∥x-ξ∥2·π(x∣u)dx.
Since the loss is strictly convex and  (11) represents a parameter integral with respect to ξ, the unique solution can be derived from
(12)∇ξ𝔼[L(X,ξ)∣U=u]=!0.
This means for every coordinate i∈{1,…,m}(13)∫ΩX‍2(xi-ξi)·π(x∣u)dx=!0,ξi·∫ΩX‍π(x∣u)dx︸=1=∫ΩX‍xi·π(x∣u)dx.
This leads to the corresponding MMSE for each coordinate i∈{1,…,m} with
(14)δi(u)=𝔼[Xi∣U=u]=c·∫ΩX‍xi·ℒu(x)·π(x)dx,
which is the expectation value of the posterior for every coordinate i∈{1,…,m}, with the normalization constant c. This example already points out the importance of the shape of the posterior π(x∣u) for the Bayesian approach.


## 2.2.4. Maximum A Posteriori (MAP)

An alternative approach, conceptually in between ML and BE, is the Maximum A Posteriori (MAP) estimation which already utilizes prior knowledge aboutx∈ΩX, but has many characteristics of the ML approach.Definition 7.
The Maximum A Posteriori estimator (MAP) for a given observationu∈ΩU is
(15)δ(u)≔argmaxxπ(x∣u)=argmaxxℒu(x)·π(x),
since the marginal distribution is positive and independent of x. In other words, MAP is the mode of the posterior.When applying the logarithm toℒu(x)·π(x), this approach can be identified with a specific Tikhonov regularization of the ML approach (also known as penalized Maximum Likelihood), as discussed in Eldar [14]. Since the MAP estimation depends on the shape of the likelihood and prior function which is problem specifically defined, under fairly general assumptions, the existence but not uniqueness of the estimate can be guaranteed. The asymptotic properties are the same as for the ML case, with some further regularity conditions on ℒu(x) and π(x); for example, see Robert [13, p. 166]. Although MAP utilizes prior knowledge about x∈ΩX, the numerical problem is the same as for ML, that is, resulting in a numerical optimization problem with the problem of identifying local and global maxima.In the context of Bayesian estimation, the MAP estimation can be derived as a limit of Bayesian estimators with respect to the constant loss function  (5). It is shown that
(16)δ(u)≔argminξ𝔼[L(X,ξ)∣U=u]=argminξ∫ΩX‍L(x,ξ)·π(x∣u)dx=argmaxξ∫Cξ‍π(x∣u)dx,
with Cξ as the m-dimensional cube with side length 2ε and midpoint ξ. In other words, the aim is to find the midpoint of the m-dimensional cube Cξ, so that the volume under the graph of the posterior is maximized. By assuming an at least continuous posterior and letting ε→0, the best midpoint δ(u) is the mode of the posterior, that is, the MAP estimator. This is illustrated in  Figure 1. Since the constant loss function is not strictly convex, the corresponding BEs are not necessarily unique, and this implicates the nonuniqueness of the MAP as a limit of these BEs.Figure 1
The MAP estimator as a limit of Bayesian estimators withx∈ℝ2. For every fixed ε the Bayesian estimator is the midpoint δ(u) which maximizes the gray shaded volume.


## 2.3. Characterization of the Estimators

The presented estimators are closely related which is graphically illustrated in  Figure2.  A major distinction is drawn in between the use of prior knowledge about x∈ΩX or not, which is shown as a thick boundary in  Figure 2. In detail, the main connections are summarized in the following list.  (i)
BE→ MMSE: the general BE reduces to the MMSE when applying the quadratic loss of  (4).(ii)
MMSE→ MAP: if the posterior distribution is symmetric, that is, π(μ+x∣u)=π(μ-x∣u) around the global maximum μ, then the MMSE estimator is identical with the MAP estimator. In this case, due to symmetry the expectation and the maximum of the posterior coincide.(iii)
BE→ MAP: the BE converges to the MAP estimator when applying the constant loss of  (5) with ε→0.(iv)
MAP→ ML: if the prior is noninformative, that is, constant for a finite ΩX or asymptotically noninformative for an infinite ΩX as the deviation parameters of the prior go to infinity, then MAP coincides with the ML estimate.Figure 2
Illustration of the relation between estimation based on ML, MAP, MMSE, and Bayesian estimation. The intersections state under which conditions the estimators are equivalent. Every condition appears only once; that is, since the conditions are not excluding each other, multiple conditions lead to further equalities.Of course, the combination of the stated properties is possible, for example, for a symmetric posterior distribution around a global maximum and a noninformative prior the ML and MMSE coincide. There are several other concepts utilized in the literature for deriving non-Bayesian estimators, that is, without (possibly noninformative) prior knowledge aboutx∈ΩX, such as the following. (i)
Unbiasedness, that is, the search of an unbiased estimator which minimizes the frequentist risk locally or uniformly. This property has been regarded as very important from the beginning of estimation theory (e.g., see Cramér [15]), although it could be shown that an unbiased estimator does not always exist (see, e.g., Lehmann and Casella [12, p. 83]), and lower risk function values (6) can be achieved when introducing a suitable bias (Eldar [16]). In this context, the Minimum Variance Unbiased Estimator (MVUE) is of special interest, since it represents the efficient unbiased estimator.(ii)
Minimaxity, that is, the minimization of the maximum value of the risk function of  (6), which typically leads to rather conservative estimates, since they can be approximated by Bayesian estimators with a prior distributions concentrated at the worst case (Robert [13, p. 67]). Unfortunately, these estimates do not always exist (Robert [13, p. 69]).(iii)
Equivariance, that is, the invariance of the estimation problem under certain transformations of the estimator. This property is important for cases in which a given estimator is applied to modified estimation situations, but the subsequent investigation of the optimality of the modified estimator should be avoided. This is important for practical cases, in which the estimator is modified routinely or a certain symmetry of the estimator is demanded by the nature of the parameters one wants to estimate. For details, see Lehmann and Casella [12, p. 147–225].


## 3. Application to Stochastic Localization in Euclidean Space

The general Bayesian framework described in the previous section will be applied to a localization problem in which 2D (3D) objects are projected onto 1D (2D) planes, and the projection signals (images) are detected for multiple views around these objects. This localization problem will be referred to as multiple view triangulation close to interfaces with a 2D or 3D observation model. In principle, multiple projections can be measured for views adjacent to each other and forming a contiguous arc or could be separated by larger angles (e.g., 90 deg). The formalism works for all sets of projections and all input data containing information about the projected coordinates of the objects. Two cases are presented in the following: linear parallel projections of 2D objects and nonlinear cone beam projection of 3D objects in a convex solution spaceΩX. All transformations are purposefully presented in Euclidean space (instead of the projective space with homogeneous coordinates [1]) in order to allow intuitive interpretations of the model and the results, as well as demonstrate direct solutions for the nonlinear 3D case, and utilize straightforwardly interpretable priors. Further, we exactly know the geometry of image acquisition (i.e., use calibrated cameras) in order to specifically focus on the noisy measurements and the incorporation of prior knowledge, instead of the challenging problem of camera calibration [1, pp. 434–457].The 2D geometry is presented in  Figure3 and shows a region which contains a feature within a larger object with interfaces, whose coordinates x∈ΩX we want to estimate. It is known that x might be at a different position, but it cannot be outside the region (e.g., is constrained by the surface of the object or regional interfaces). Furthermore, the PDFs are assumed to be Gaussians or truncated Gaussians (or uniform distributions in the Results section), since these are valid assumptions in many applications [1], and the resulting estimators can partly be presented in terms of explicit relationships. In general, the estimation we present works for arbitrary PDFs and can be directly applied to nonlinear/non-Gaussian problems.Figure 3
Geometrical projectionu∈ℝ of the coordinate x∈ℝ2 in the detector plane.Constraints due to interfaces or surfaces, for example, might occur in radiotherapy (or interventional radiology, image-guided surgical procedure) where the location of an implanted fiducial marker or another feature has to be determined in real-time based on X-ray projections (radiographs or tomosynthesis projections) [3, 4]. Let us focus on a specific application of breast cancer patient setup in radiotherapy. In this case, surgical clips implanted during tumor resection are used for localization of the tumor bed. Due to deformations of the breast in between the original CT images and the treatment, the markers might be at arbitrary locations, but they will never be outside of the patient's body, and they will not be inside the chest defined by the ribs. Also it is likely that they will not be displaced by more than a given range known from clinical practice (e.g., 2 cm from the original location in the planning CT scan). In radiotherapy, the goal is to localize the tumor (or the tumor surrogate such as the marker) with 1 mm accuracy. This is achievable by performing a CT scan. But such a CT scan is not possible in many clinical situations either due to undesired dose from the imaging or due to time constraints or other practical concerns. Additionally, in the case of a given marker within the flexible breast tissue, the surface can be measured with optical or infrared cameras in real-time or based on a tomosynthesis reconstruction, which is an intensely discussed topic [17–19]. We can thus assume that the knowledge about the breast surface is available for the given posture of the patient. In a similar context, in interventional radiology or image based surgeries there is a need to work with a few projections or even single images in order to find the location of certain objects in 3D. Many other applications are thinkable, such as robotic arm position estimation in a closed room or computer vision approaches with spatial restrictions [1].


### 3.1. Parallel Beam in 2D

In this casex∈ΩX⊂ℝ2. The observation u∈ℝ is performed by projecting the feature in a detector plane at different angles Θ. The aim of this example is to estimate the coordinate x based on the observed coordinates u at N different angles between 0∘ and 90∘, which are merged to a single vector u=(u1,…,uN)T∈ΩU⊂ℝN.Please note that although the presented situation might have several applications, the numerical example in the Results section uses arbitrary dimensionless parameters in order to demonstrate the principal differences of the estimation methods. Furthermore, the linear relation between observation and the coordinate in the parallel beam model allows to derive explicit estimators, which is not possible for the nonlinear cone beam geometry in  Section3.2.


#### 3.1.1. Observation without Uncertainties and Prior Knowledge

If it is assumed that the observation is without any errors, then no prior knowledge is necessary because there is no ambiguity in the detection, and the linear relation betweenx and a single observation u can be derived from
(17)ui≔-x1sinΘi+x2cosΘi.
This is a linear equation in the two coordinates of x and; therefore, two equations need to be derived for a unique solution, that is, two observations u1 and u2 at angles Θ1 and Θ2 with Θ1≠Θ2 modulo 360∘. The explicit solution for this straightforward estimation is given by the estimator δ(u)(18)δ(u)≔(-sinΘ1cosΘ1-sinΘ2cosΘ2)-1·(u1u2).
In the numerical examples we will assume the two angles to be Θ1=0∘ and Θ2=90∘ which represents the best separation between two projections. In the further text the estimator of  (18) will be denoted by the 2×2 approach.


#### 3.1.2. Observation with Uncertainties and without Prior Knowledge

If we know that there are observational uncertainties, that is, errors in determiningu, then one can assume an error distribution, which is represented in this example by the normally distributed random variable at angle Θi(19)ΔUi~𝒩(0,σΔU2).
In the geometrical model, this means we are adding a random error to the observed projected coordinate at every angle
(20)Ui≔-x1sinΘi+x2cosΘi+ΔUi,
which makes now the observation itself to a random variable with the PDF
(21)f[Ui](ui∣x)=1σΔU2π×exp[-(ui-(-x1sinΘi+x2cosΘi))22σΔU2].
By acquiring N≥2 stochastically independent observations u at angles Θ1,…,ΘN, one can summarize these observations in a single PDF, which is the 2D likelihood function in x(22)f[U](u∣x)=∏i=1N‍f[Ui](ui∣x)=1(σΔU2π)N×exp[-∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))22σΔU2].
In consequence, by applying the logarithm the Maximum Likelihood estimator is determined by
(23)δ(u)≔argminx∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))2,
which is exactly the least squares solution. This is a well-known result in estimation theory and indicates the success of least squares, since the sum of many independent errors can often be approximated well by a normal distribution [1, pp. 102, 121, 134, 434, and 568]. An important observation is that in this approach the solution is totally independent of the knowledge about the measurement uncertainty σΔU. The explicit solution is derived from solving (23) analytically, and it follows (cp.  (18) )
(24)δ(u)=(∑i=1Nsin2Θi-∑i=1N‍sinΘicosΘi-∑i=1N‍sinΘicosΘi∑i=1N‍cos2Θi)-1·(-∑i=1N‍uisinΘi∑i=1N‍uicosΘi).


#### 3.1.3. Observation with Uncertainties and Prior Knowledge

If it is further known thatx is stochastically distributed, that is, a random variable X, one can additionally introduce prior knowledge about this distribution into the estimation model. In the 2D example the prior is defined as follows
(25)X~π(x)≔a·𝒩(μX,ΣX)(x)·χC(m,r)(x),
with a normalization constant a, the characteristic function χ which has the value 1 in the circle C(m,r) with midpoint m and radius r, and 0 elsewhere (representing an example of a convex region), and μX∈ℝ2 the mean of the normal distribution, and
(26)ΣX≔(σX1200σX22).
This is essentially a Gaussian density function which is truncated by the circle C(m,r), which is illustrated in  Figure 4.  In the geometrical model now also the parameters become random variables; that is,
(27)Ui≔-X1sinΘi+X2cosΘi+ΔUi.
Utilizing Bayes’ theorem, the posterior PDF can be calculated by (c is the normalization constant of the posterior)
(28)π(x∣u)=c·f[U](u∣x)·π(x)=c·∏i=1N‍f[Ui](ui∣x)·π(x)(29)π(x∣u)=c~·exp[(xj-μXj)22σXj2-∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))22σΔU2==c~·exp-∑j=12‍(xj-μXj)22σXj2]·χC(m,r)(x),(30)c~≔ca(σΔU)NσX1σX2(2π)N/2+1>0.
Based on the posterior the Maximum A Posteriori estimator is therefore defined as
(31)δ(u)≔argminx∈C(m,r)∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))2σΔU2+∑j=12‍(xj-μXj)2σXj2.
This formula essentially shows explicitly how the compromise between measurement and prior knowledge, as two different sources of information, is gained: the measurement is weighted with 1/σΔU2 against the prior knowledge with 1/σXj2 (j=1,2), which is the case of the Tikhonov regularization for the constrained least squares problem. Since this is a convex function on a convex domain C(m,r), there is a guaranteed global minimum. If this maximum is in the interior of C(m,r), the explicit solution is determined by solving  (31), and it follows (cp.  (18) and  (24))
(32)δ(u)=(∑i=1N‍sin2Θi+(σΔUσX1)2-∑i=1N‍sinΘicosΘi-∑i=1N‍sinΘicosΘi∑i=1N‍cos2Θi+(σΔUσX2)2)-1·(-∑i=1N‍uisinΘi+μX1(σΔUσX1)2∑i=1N‍uicosΘi+μX2(σΔUσX1)2),
From a practical perspective, if the δ(u) calculated by  (32) is outside of C(m,r), then the maximum of  (31) must be on the border (circle in  Figure 4), which can then be solved by Lagrangian multipliers applied to  (31) or directly by solving  (31) with a general optimization algorithm. The explicit solution for the interior of C(m,r) depends now on the ratios of σΔU/σXi (i∈{1,2}), which again shows a weighting of the measurement uncertainties and the prior knowledge. Particularly, if σX1/2→∞, the prior gets noninformative, and the interior estimator converges to the ML estimator.(a) Example for the prior density function with the parameters specified in the Results section. (b) 10,000 samples drawn from the prior distribution which are utilized in the Results section.
(a)
(b)Based on the posterior also the MinimumMean Square Error estimator can be determined by
(33)δi(u)≔𝔼[Xi∣u]=∬C(m,r)xi·π(x∣u)dx
for i∈{1,2} with the posterior in  (29), which cannot be calculated elementarily anymore. For the specific case of r→∞, the MMSE estimator is coinciding with the MAP estimator of  (32) due to the symmetry of the posterior around its global maximum when using Gaussians. The difference between MAP and MMSE will therefore be enhanced by the estimation close to the boundary. Please note that both MAP and MMSE guarantee an estimate which is not outside of the (convex) circle C(m,r), effectively bounding the feasibility region with the prior. In contrast, for the 2×2 and ML approach no bounded feasibility region is possible.


### 3.2. Cone Beam in 3D

In this casex∈ΩX⊂ℝ3. The observation u∈ℝ2 is again performed by projecting the feature in a detector plane at different angles Θ. There are N observed coordinates u at different angles between 0∘ and 90∘, which are merged to a single matrix u=(u1,…,uN)T∈ΩU⊂ℝ2×N (which can be identified with a vector in ℝ2N).If it is assumed that the observation is without error and no prior knowledge is available to the estimation, the nonlinear relation betweenx and a single observation ui at angle Θi is
(34)ui=ψi(x)≔g+hx1cosΘi+x2sinΘi+g×(-x1sinΘi+x2cosΘix3),
with the imaging source to isocenter of rotation distance g and isocenter to detection plane distance h; see  Figure 5 for an illustration of this geometry. This projection excludes points x which are on the plane x1cosΘi+x2sinΘi+g=0 orthogonal to the beam through the imaging source. It is directly recognizable that ψi is nonlinear (in Euclidean space), since ψi(x+y)≠ψi(x)+ψi(y). We recognize that this projection can be represented as a linear transformation in homogeneous coordinates (e.g., see Hartley and Zisserman [1]), but we purposefully ignore this in order to show how one can apply the proposed estimation approaches directly to a nonlinear situation as well as define intuitive priors in Euclidean space.Figure 5
Geometrical projectionu∈ℝ2 of the coordinate x∈ℝ3 in the detector plane.A direct solution, such as the2×2 case in 2D, is not uniquely possible since there are equation numbers as a multiple of 2, but we have 3 unknowns. We avoid this ambiguity by not defining a corresponding 3×3 case.


#### 3.2.1. Observation with Uncertainties and without Prior Knowledge

The observation uncertainty at angleΘi is in this case described by
(35)ΔUi~𝒩(0,σΔU·I2×2),
with the 2×2 identity matrix I2×2. The geometrical model is therefore modified to
(36)Ui≔ψi(x)+ΔUi.
In consequence, the observation PDF is
(37)f[Ui](ui∣x)=1σΔU22πexp[-∑k=12‍(uk,i-ψk,i(x))22σΔU2].
By acquiring N≥2 stochastically independent observations ui at angles Θ1,…,ΘN, one can summarize these observation in a single PDF, which is the 3D likelihood function in x(38)f[U](u∣x)=1(σΔU22π)N×exp[-∑i=1N‍∑k=12‍(uk,i-ψk,i(x))22σΔU2].
This leads then to the Maximum Likelihood estimator
(39)δ(u)≔argminx×∑i=1N‍(u1,i-(-x1sinΘi+x2cosΘi)(g+h)x1cosΘi+x2sinΘi+g)2+(u2,i-x3(g+h)x1cosΘi+x2sinΘi+g)2,
which is the nonlinear least squares solution. Due to the nonlinearity, no elementary solution can be derived, and this expression is directly minimized.


#### 3.2.2. Observation with Uncertainties and Prior Knowledge

The prior knowledge is in the 3D case modeled with(40)X~π(x)≔a·𝒩(μX,ΣX)(x)·χS(m,r)(x),
with a normalization constant a, the characteristic function χ which has the value 1 in the sphere S(m,r) with midpoint m and radius r, and 0 elsewhere (again as an example of a convex region), and μX∈ℝ3 the mean of the normal distribution, and
(41)ΣX≔(σX12000σX22000σX32).
This is essentially a Gaussian density function which is truncated by the sphere S(m,r). A sample from this prior is illustrated in  Figure 6. The geometrical model is then
(42)Ui≔ψi(X)+ΔUi.
The posterior distribution then gets (c is the normalization constant of the posterior)
(43)π(x∣u)=c~·exp[-∑i=1N‍∑k=12‍(uk,i-ψk,i(x))22σΔU2-∑j=13‍(xj-μXj)22σXj2]·χS(m,r)(x),(44)c~≔ca(σΔU)2NσX1σX2σX3(2π)N+3/2>0.
In consequence, the Maximum A Posteriori estimator is
(45)δ(u)≔argminx∈S(m,r)×∑i=1N‍((u1,i-(-x1sinΘi+x2cosΘi)(g+h)x1cosΘi+x2sinΘi+g)2+(u2,i-x3(g+h)x1cosΘi+x2sinΘi+g)2)×(σΔU2)-1+∑j=13‍(xj-μXj)2σXj2,
representing a constrained and regularized nonlinear least squares optimization problem. The corresponding Minimum Mean Square Error estimator is analogously
(46)δi(u)≔𝔼[Xi∣u]=∭S(m,r)xi·π(x∣u)dx,
for i∈{1,2,3} with the posterior of  (43), which again cannot be calculated elementarily anymore.Figure 6
10,000 samples drawn from the prior distribution which are utilized in the Results section.


## 3.1. Parallel Beam in 2D

In this casex∈ΩX⊂ℝ2. The observation u∈ℝ is performed by projecting the feature in a detector plane at different angles Θ. The aim of this example is to estimate the coordinate x based on the observed coordinates u at N different angles between 0∘ and 90∘, which are merged to a single vector u=(u1,…,uN)T∈ΩU⊂ℝN.Please note that although the presented situation might have several applications, the numerical example in the Results section uses arbitrary dimensionless parameters in order to demonstrate the principal differences of the estimation methods. Furthermore, the linear relation between observation and the coordinate in the parallel beam model allows to derive explicit estimators, which is not possible for the nonlinear cone beam geometry in  Section3.2.


### 3.1.1. Observation without Uncertainties and Prior Knowledge

If it is assumed that the observation is without any errors, then no prior knowledge is necessary because there is no ambiguity in the detection, and the linear relation betweenx and a single observation u can be derived from
(17)ui≔-x1sinΘi+x2cosΘi.
This is a linear equation in the two coordinates of x and; therefore, two equations need to be derived for a unique solution, that is, two observations u1 and u2 at angles Θ1 and Θ2 with Θ1≠Θ2 modulo 360∘. The explicit solution for this straightforward estimation is given by the estimator δ(u)(18)δ(u)≔(-sinΘ1cosΘ1-sinΘ2cosΘ2)-1·(u1u2).
In the numerical examples we will assume the two angles to be Θ1=0∘ and Θ2=90∘ which represents the best separation between two projections. In the further text the estimator of  (18) will be denoted by the 2×2 approach.


### 3.1.2. Observation with Uncertainties and without Prior Knowledge

If we know that there are observational uncertainties, that is, errors in determiningu, then one can assume an error distribution, which is represented in this example by the normally distributed random variable at angle Θi(19)ΔUi~𝒩(0,σΔU2).
In the geometrical model, this means we are adding a random error to the observed projected coordinate at every angle
(20)Ui≔-x1sinΘi+x2cosΘi+ΔUi,
which makes now the observation itself to a random variable with the PDF
(21)f[Ui](ui∣x)=1σΔU2π×exp[-(ui-(-x1sinΘi+x2cosΘi))22σΔU2].
By acquiring N≥2 stochastically independent observations u at angles Θ1,…,ΘN, one can summarize these observations in a single PDF, which is the 2D likelihood function in x(22)f[U](u∣x)=∏i=1N‍f[Ui](ui∣x)=1(σΔU2π)N×exp[-∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))22σΔU2].
In consequence, by applying the logarithm the Maximum Likelihood estimator is determined by
(23)δ(u)≔argminx∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))2,
which is exactly the least squares solution. This is a well-known result in estimation theory and indicates the success of least squares, since the sum of many independent errors can often be approximated well by a normal distribution [1, pp. 102, 121, 134, 434, and 568]. An important observation is that in this approach the solution is totally independent of the knowledge about the measurement uncertainty σΔU. The explicit solution is derived from solving (23) analytically, and it follows (cp.  (18) )
(24)δ(u)=(∑i=1Nsin2Θi-∑i=1N‍sinΘicosΘi-∑i=1N‍sinΘicosΘi∑i=1N‍cos2Θi)-1·(-∑i=1N‍uisinΘi∑i=1N‍uicosΘi).


### 3.1.3. Observation with Uncertainties and Prior Knowledge

If it is further known thatx is stochastically distributed, that is, a random variable X, one can additionally introduce prior knowledge about this distribution into the estimation model. In the 2D example the prior is defined as follows
(25)X~π(x)≔a·𝒩(μX,ΣX)(x)·χC(m,r)(x),
with a normalization constant a, the characteristic function χ which has the value 1 in the circle C(m,r) with midpoint m and radius r, and 0 elsewhere (representing an example of a convex region), and μX∈ℝ2 the mean of the normal distribution, and
(26)ΣX≔(σX1200σX22).
This is essentially a Gaussian density function which is truncated by the circle C(m,r), which is illustrated in  Figure 4.  In the geometrical model now also the parameters become random variables; that is,
(27)Ui≔-X1sinΘi+X2cosΘi+ΔUi.
Utilizing Bayes’ theorem, the posterior PDF can be calculated by (c is the normalization constant of the posterior)
(28)π(x∣u)=c·f[U](u∣x)·π(x)=c·∏i=1N‍f[Ui](ui∣x)·π(x)(29)π(x∣u)=c~·exp[(xj-μXj)22σXj2-∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))22σΔU2==c~·exp-∑j=12‍(xj-μXj)22σXj2]·χC(m,r)(x),(30)c~≔ca(σΔU)NσX1σX2(2π)N/2+1>0.
Based on the posterior the Maximum A Posteriori estimator is therefore defined as
(31)δ(u)≔argminx∈C(m,r)∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))2σΔU2+∑j=12‍(xj-μXj)2σXj2.
This formula essentially shows explicitly how the compromise between measurement and prior knowledge, as two different sources of information, is gained: the measurement is weighted with 1/σΔU2 against the prior knowledge with 1/σXj2 (j=1,2), which is the case of the Tikhonov regularization for the constrained least squares problem. Since this is a convex function on a convex domain C(m,r), there is a guaranteed global minimum. If this maximum is in the interior of C(m,r), the explicit solution is determined by solving  (31), and it follows (cp.  (18) and  (24))
(32)δ(u)=(∑i=1N‍sin2Θi+(σΔUσX1)2-∑i=1N‍sinΘicosΘi-∑i=1N‍sinΘicosΘi∑i=1N‍cos2Θi+(σΔUσX2)2)-1·(-∑i=1N‍uisinΘi+μX1(σΔUσX1)2∑i=1N‍uicosΘi+μX2(σΔUσX1)2),
From a practical perspective, if the δ(u) calculated by  (32) is outside of C(m,r), then the maximum of  (31) must be on the border (circle in  Figure 4), which can then be solved by Lagrangian multipliers applied to  (31) or directly by solving  (31) with a general optimization algorithm. The explicit solution for the interior of C(m,r) depends now on the ratios of σΔU/σXi (i∈{1,2}), which again shows a weighting of the measurement uncertainties and the prior knowledge. Particularly, if σX1/2→∞, the prior gets noninformative, and the interior estimator converges to the ML estimator.(a) Example for the prior density function with the parameters specified in the Results section. (b) 10,000 samples drawn from the prior distribution which are utilized in the Results section.
(a)
(b)Based on the posterior also the MinimumMean Square Error estimator can be determined by
(33)δi(u)≔𝔼[Xi∣u]=∬C(m,r)xi·π(x∣u)dx
for i∈{1,2} with the posterior in  (29), which cannot be calculated elementarily anymore. For the specific case of r→∞, the MMSE estimator is coinciding with the MAP estimator of  (32) due to the symmetry of the posterior around its global maximum when using Gaussians. The difference between MAP and MMSE will therefore be enhanced by the estimation close to the boundary. Please note that both MAP and MMSE guarantee an estimate which is not outside of the (convex) circle C(m,r), effectively bounding the feasibility region with the prior. In contrast, for the 2×2 and ML approach no bounded feasibility region is possible.


## 3.1.1. Observation without Uncertainties and Prior Knowledge

If it is assumed that the observation is without any errors, then no prior knowledge is necessary because there is no ambiguity in the detection, and the linear relation betweenx and a single observation u can be derived from
(17)ui≔-x1sinΘi+x2cosΘi.
This is a linear equation in the two coordinates of x and; therefore, two equations need to be derived for a unique solution, that is, two observations u1 and u2 at angles Θ1 and Θ2 with Θ1≠Θ2 modulo 360∘. The explicit solution for this straightforward estimation is given by the estimator δ(u)(18)δ(u)≔(-sinΘ1cosΘ1-sinΘ2cosΘ2)-1·(u1u2).
In the numerical examples we will assume the two angles to be Θ1=0∘ and Θ2=90∘ which represents the best separation between two projections. In the further text the estimator of  (18) will be denoted by the 2×2 approach.


## 3.1.2. Observation with Uncertainties and without Prior Knowledge

If we know that there are observational uncertainties, that is, errors in determiningu, then one can assume an error distribution, which is represented in this example by the normally distributed random variable at angle Θi(19)ΔUi~𝒩(0,σΔU2).
In the geometrical model, this means we are adding a random error to the observed projected coordinate at every angle
(20)Ui≔-x1sinΘi+x2cosΘi+ΔUi,
which makes now the observation itself to a random variable with the PDF
(21)f[Ui](ui∣x)=1σΔU2π×exp[-(ui-(-x1sinΘi+x2cosΘi))22σΔU2].
By acquiring N≥2 stochastically independent observations u at angles Θ1,…,ΘN, one can summarize these observations in a single PDF, which is the 2D likelihood function in x(22)f[U](u∣x)=∏i=1N‍f[Ui](ui∣x)=1(σΔU2π)N×exp[-∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))22σΔU2].
In consequence, by applying the logarithm the Maximum Likelihood estimator is determined by
(23)δ(u)≔argminx∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))2,
which is exactly the least squares solution. This is a well-known result in estimation theory and indicates the success of least squares, since the sum of many independent errors can often be approximated well by a normal distribution [1, pp. 102, 121, 134, 434, and 568]. An important observation is that in this approach the solution is totally independent of the knowledge about the measurement uncertainty σΔU. The explicit solution is derived from solving (23) analytically, and it follows (cp.  (18) )
(24)δ(u)=(∑i=1Nsin2Θi-∑i=1N‍sinΘicosΘi-∑i=1N‍sinΘicosΘi∑i=1N‍cos2Θi)-1·(-∑i=1N‍uisinΘi∑i=1N‍uicosΘi).


## 3.1.3. Observation with Uncertainties and Prior Knowledge

If it is further known thatx is stochastically distributed, that is, a random variable X, one can additionally introduce prior knowledge about this distribution into the estimation model. In the 2D example the prior is defined as follows
(25)X~π(x)≔a·𝒩(μX,ΣX)(x)·χC(m,r)(x),
with a normalization constant a, the characteristic function χ which has the value 1 in the circle C(m,r) with midpoint m and radius r, and 0 elsewhere (representing an example of a convex region), and μX∈ℝ2 the mean of the normal distribution, and
(26)ΣX≔(σX1200σX22).
This is essentially a Gaussian density function which is truncated by the circle C(m,r), which is illustrated in  Figure 4.  In the geometrical model now also the parameters become random variables; that is,
(27)Ui≔-X1sinΘi+X2cosΘi+ΔUi.
Utilizing Bayes’ theorem, the posterior PDF can be calculated by (c is the normalization constant of the posterior)
(28)π(x∣u)=c·f[U](u∣x)·π(x)=c·∏i=1N‍f[Ui](ui∣x)·π(x)(29)π(x∣u)=c~·exp[(xj-μXj)22σXj2-∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))22σΔU2==c~·exp-∑j=12‍(xj-μXj)22σXj2]·χC(m,r)(x),(30)c~≔ca(σΔU)NσX1σX2(2π)N/2+1>0.
Based on the posterior the Maximum A Posteriori estimator is therefore defined as
(31)δ(u)≔argminx∈C(m,r)∑i=1N‍(ui-(-x1sinΘi+x2cosΘi))2σΔU2+∑j=12‍(xj-μXj)2σXj2.
This formula essentially shows explicitly how the compromise between measurement and prior knowledge, as two different sources of information, is gained: the measurement is weighted with 1/σΔU2 against the prior knowledge with 1/σXj2 (j=1,2), which is the case of the Tikhonov regularization for the constrained least squares problem. Since this is a convex function on a convex domain C(m,r), there is a guaranteed global minimum. If this maximum is in the interior of C(m,r), the explicit solution is determined by solving  (31), and it follows (cp.  (18) and  (24))
(32)δ(u)=(∑i=1N‍sin2Θi+(σΔUσX1)2-∑i=1N‍sinΘicosΘi-∑i=1N‍sinΘicosΘi∑i=1N‍cos2Θi+(σΔUσX2)2)-1·(-∑i=1N‍uisinΘi+μX1(σΔUσX1)2∑i=1N‍uicosΘi+μX2(σΔUσX1)2),
From a practical perspective, if the δ(u) calculated by  (32) is outside of C(m,r), then the maximum of  (31) must be on the border (circle in  Figure 4), which can then be solved by Lagrangian multipliers applied to  (31) or directly by solving  (31) with a general optimization algorithm. The explicit solution for the interior of C(m,r) depends now on the ratios of σΔU/σXi (i∈{1,2}), which again shows a weighting of the measurement uncertainties and the prior knowledge. Particularly, if σX1/2→∞, the prior gets noninformative, and the interior estimator converges to the ML estimator.(a) Example for the prior density function with the parameters specified in the Results section. (b) 10,000 samples drawn from the prior distribution which are utilized in the Results section.
(a)
(b)Based on the posterior also the MinimumMean Square Error estimator can be determined by
(33)δi(u)≔𝔼[Xi∣u]=∬C(m,r)xi·π(x∣u)dx
for i∈{1,2} with the posterior in  (29), which cannot be calculated elementarily anymore. For the specific case of r→∞, the MMSE estimator is coinciding with the MAP estimator of  (32) due to the symmetry of the posterior around its global maximum when using Gaussians. The difference between MAP and MMSE will therefore be enhanced by the estimation close to the boundary. Please note that both MAP and MMSE guarantee an estimate which is not outside of the (convex) circle C(m,r), effectively bounding the feasibility region with the prior. In contrast, for the 2×2 and ML approach no bounded feasibility region is possible.


## 3.2. Cone Beam in 3D

In this casex∈ΩX⊂ℝ3. The observation u∈ℝ2 is again performed by projecting the feature in a detector plane at different angles Θ. There are N observed coordinates u at different angles between 0∘ and 90∘, which are merged to a single matrix u=(u1,…,uN)T∈ΩU⊂ℝ2×N (which can be identified with a vector in ℝ2N).If it is assumed that the observation is without error and no prior knowledge is available to the estimation, the nonlinear relation betweenx and a single observation ui at angle Θi is
(34)ui=ψi(x)≔g+hx1cosΘi+x2sinΘi+g×(-x1sinΘi+x2cosΘix3),
with the imaging source to isocenter of rotation distance g and isocenter to detection plane distance h; see  Figure 5 for an illustration of this geometry. This projection excludes points x which are on the plane x1cosΘi+x2sinΘi+g=0 orthogonal to the beam through the imaging source. It is directly recognizable that ψi is nonlinear (in Euclidean space), since ψi(x+y)≠ψi(x)+ψi(y). We recognize that this projection can be represented as a linear transformation in homogeneous coordinates (e.g., see Hartley and Zisserman [1]), but we purposefully ignore this in order to show how one can apply the proposed estimation approaches directly to a nonlinear situation as well as define intuitive priors in Euclidean space.Figure 5
Geometrical projectionu∈ℝ2 of the coordinate x∈ℝ3 in the detector plane.A direct solution, such as the2×2 case in 2D, is not uniquely possible since there are equation numbers as a multiple of 2, but we have 3 unknowns. We avoid this ambiguity by not defining a corresponding 3×3 case.


### 3.2.1. Observation with Uncertainties and without Prior Knowledge

The observation uncertainty at angleΘi is in this case described by
(35)ΔUi~𝒩(0,σΔU·I2×2),
with the 2×2 identity matrix I2×2. The geometrical model is therefore modified to
(36)Ui≔ψi(x)+ΔUi.
In consequence, the observation PDF is
(37)f[Ui](ui∣x)=1σΔU22πexp[-∑k=12‍(uk,i-ψk,i(x))22σΔU2].
By acquiring N≥2 stochastically independent observations ui at angles Θ1,…,ΘN, one can summarize these observation in a single PDF, which is the 3D likelihood function in x(38)f[U](u∣x)=1(σΔU22π)N×exp[-∑i=1N‍∑k=12‍(uk,i-ψk,i(x))22σΔU2].
This leads then to the Maximum Likelihood estimator
(39)δ(u)≔argminx×∑i=1N‍(u1,i-(-x1sinΘi+x2cosΘi)(g+h)x1cosΘi+x2sinΘi+g)2+(u2,i-x3(g+h)x1cosΘi+x2sinΘi+g)2,
which is the nonlinear least squares solution. Due to the nonlinearity, no elementary solution can be derived, and this expression is directly minimized.


### 3.2.2. Observation with Uncertainties and Prior Knowledge

The prior knowledge is in the 3D case modeled with(40)X~π(x)≔a·𝒩(μX,ΣX)(x)·χS(m,r)(x),
with a normalization constant a, the characteristic function χ which has the value 1 in the sphere S(m,r) with midpoint m and radius r, and 0 elsewhere (again as an example of a convex region), and μX∈ℝ3 the mean of the normal distribution, and
(41)ΣX≔(σX12000σX22000σX32).
This is essentially a Gaussian density function which is truncated by the sphere S(m,r). A sample from this prior is illustrated in  Figure 6. The geometrical model is then
(42)Ui≔ψi(X)+ΔUi.
The posterior distribution then gets (c is the normalization constant of the posterior)
(43)π(x∣u)=c~·exp[-∑i=1N‍∑k=12‍(uk,i-ψk,i(x))22σΔU2-∑j=13‍(xj-μXj)22σXj2]·χS(m,r)(x),(44)c~≔ca(σΔU)2NσX1σX2σX3(2π)N+3/2>0.
In consequence, the Maximum A Posteriori estimator is
(45)δ(u)≔argminx∈S(m,r)×∑i=1N‍((u1,i-(-x1sinΘi+x2cosΘi)(g+h)x1cosΘi+x2sinΘi+g)2+(u2,i-x3(g+h)x1cosΘi+x2sinΘi+g)2)×(σΔU2)-1+∑j=13‍(xj-μXj)2σXj2,
representing a constrained and regularized nonlinear least squares optimization problem. The corresponding Minimum Mean Square Error estimator is analogously
(46)δi(u)≔𝔼[Xi∣u]=∭S(m,r)xi·π(x∣u)dx,
for i∈{1,2,3} with the posterior of  (43), which again cannot be calculated elementarily anymore.Figure 6
10,000 samples drawn from the prior distribution which are utilized in the Results section.


## 3.2.1. Observation with Uncertainties and without Prior Knowledge

The observation uncertainty at angleΘi is in this case described by
(35)ΔUi~𝒩(0,σΔU·I2×2),
with the 2×2 identity matrix I2×2. The geometrical model is therefore modified to
(36)Ui≔ψi(x)+ΔUi.
In consequence, the observation PDF is
(37)f[Ui](ui∣x)=1σΔU22πexp[-∑k=12‍(uk,i-ψk,i(x))22σΔU2].
By acquiring N≥2 stochastically independent observations ui at angles Θ1,…,ΘN, one can summarize these observation in a single PDF, which is the 3D likelihood function in x(38)f[U](u∣x)=1(σΔU22π)N×exp[-∑i=1N‍∑k=12‍(uk,i-ψk,i(x))22σΔU2].
This leads then to the Maximum Likelihood estimator
(39)δ(u)≔argminx×∑i=1N‍(u1,i-(-x1sinΘi+x2cosΘi)(g+h)x1cosΘi+x2sinΘi+g)2+(u2,i-x3(g+h)x1cosΘi+x2sinΘi+g)2,
which is the nonlinear least squares solution. Due to the nonlinearity, no elementary solution can be derived, and this expression is directly minimized.


## 3.2.2. Observation with Uncertainties and Prior Knowledge

The prior knowledge is in the 3D case modeled with(40)X~π(x)≔a·𝒩(μX,ΣX)(x)·χS(m,r)(x),
with a normalization constant a, the characteristic function χ which has the value 1 in the sphere S(m,r) with midpoint m and radius r, and 0 elsewhere (again as an example of a convex region), and μX∈ℝ3 the mean of the normal distribution, and
(41)ΣX≔(σX12000σX22000σX32).
This is essentially a Gaussian density function which is truncated by the sphere S(m,r). A sample from this prior is illustrated in  Figure 6. The geometrical model is then
(42)Ui≔ψi(X)+ΔUi.
The posterior distribution then gets (c is the normalization constant of the posterior)
(43)π(x∣u)=c~·exp[-∑i=1N‍∑k=12‍(uk,i-ψk,i(x))22σΔU2-∑j=13‍(xj-μXj)22σXj2]·χS(m,r)(x),(44)c~≔ca(σΔU)2NσX1σX2σX3(2π)N+3/2>0.
In consequence, the Maximum A Posteriori estimator is
(45)δ(u)≔argminx∈S(m,r)×∑i=1N‍((u1,i-(-x1sinΘi+x2cosΘi)(g+h)x1cosΘi+x2sinΘi+g)2+(u2,i-x3(g+h)x1cosΘi+x2sinΘi+g)2)×(σΔU2)-1+∑j=13‍(xj-μXj)2σXj2,
representing a constrained and regularized nonlinear least squares optimization problem. The corresponding Minimum Mean Square Error estimator is analogously
(46)δi(u)≔𝔼[Xi∣u]=∭S(m,r)xi·π(x∣u)dx,
for i∈{1,2,3} with the posterior of  (43), which again cannot be calculated elementarily anymore.Figure 6
10,000 samples drawn from the prior distribution which are utilized in the Results section.


## 4. Results

In order to show the practical differences, we present 5 different cases for the comparison of the estimators with 10,000 sample pointsx each (recall that X is a random variable), in 2D parallel beam and 3D cone beam geometry. The standard parameter set (A) is as follows: with N=5 projections equidistantly distributed in between 0∘ and 90∘ and σΔU=3, σX1=σX2=σX3=3, μX=(16.5,16.5)T, or μX=(16.5,16.5,16.5)T, respectively, m=(10,10)T or m=(10,10,10)T respectively, r=10,  g=100 and  h=50.Besides the standard parameter set presented above, four other parameter sets (cases) are investigated (by changing only one parameter at a time while setting all other parameters to the standard): (B) applying onlyN=2 projections, (C) applying N=10 projections, (D) narrowing the prior by σXi=1.5, and (E) reducing the measurement uncertainties σΔU=1.5. In addition, for the 2D parallel beam case, we present the results of the MAP and MMSE estimations if there is a noninformative prior (np) utilized, that is, uniform distribution inside the circle C(m,r) and zero outside,
(47)σXi⟶∞,i∈{1,2}⟹π(x)≔1r2π·χC(m,r)(x),
which are denoted with MAPnp≔limσXi→∞MAP and MMSEnp accordingly. This means, only the hard constraints about the feasibility region are utilized in these estimations, without prior knowledge about the concentration of x at the border (circle in  Figure 4). In a broader context, such a noninformative MAPnp estimator can be identified with a ML estimator that has ad hoc constraints about the solution space in optimization.


### 4.1. Parallel Beam in 2D

For the standard parameters of case (A), in  Figure7 an example of the posterior PDF in 2D for a single sample point x is presented. In this plot on the right side a profile of the likelihood and the posterior function is presented, which demonstrates the general behavior: the width of the posterior is reduced compared to the likelihood function by the cost of introducing a shift (bias) of the maximum position by prior knowledge. Further, in the posterior profile the circular boundary is presented by the truncation of the profile. By incorporating the truncation into the prior distribution, the MMSE (which for simplicity one can think of as the expectation value of the posterior profile) is compensating for the bias of the maximum position in MAP.(a) Posterior density function for a specific estimation of case (A) and full prior. (b) Profile of the likelihood and posterior function atx=r·(cosπ/6,sinπ/6)T with r∈[0,12].
(a)
(b)In  Figure8 the estimation performance, for case (A) of the 2×2 approach of  (18), the ML of  (24), the MAP of  (32) and  (31), and the MMSE of  (33) (with noninformative prior and truncated Gaussian prior) are presented. The complete estimation results for all five cases (A)–(E) are presented in  Table 2.  In addition to  Table 2 the RMSEs of a broader variation of the essential parameters are illustrated in  Figure 9 by starting with case (A) and varying only the specified parameter of each subfigure. Examples for the posterior distributions  cases (A)–(E) (for the same specific sample x) are presented in  Figure 10, in order to see the qualitative differences due to the according changes in the parameter set.Table 2
Statistical values about the estimation results of 2D parallel beam geometry for cases (A)–(E) withM= 10,000 sample points x each are presented for the six 2D estimators. The radial estimation error is presented at first and the coordinate estimation errors of (x1,x2) are shown in brackets. In the following the formulas for the statistical values are presented for the radial error: the mean value=∥δ(u)∥-∥x∥¯=1/M∑j=1M‍∥δ(uj)∥-∥xj∥, the standard deviation=(1/(M-1)∑j=1M‍(∥δ(uj)∥-∥δ(u)∥¯-∥xj∥+∥x∥¯)2)1/2, the absolute maximum=maxj=1⋯M|∥δ(uj)∥-∥xj∥|, the Root Mean Square Error (RMSE) =(1/M∑j=1M‍(∥δ(uj)∥-∥xj∥2)1/2 and thegainof methodX relative to ML =(RMSEML-RMSEX)/RMSEML.
Case
Method
Mean
st. deviation
abs. maximum
Radial RMSE
[
%
]
2
×
2
3.77 (−0.0,0.0)
1.94 (3.0,3.0)
11.90 (11.1,10.9)
4.24 (3.0,3.0)
-
39
ML
2.67 (−0.0,0.0)
1.48 (2.2,2.2)
9.75 (9.2,8.3)
3.05 (2.2,2.2)
—
A
MAPnp
2.35 (−0.3,−0.3)
1.35 (1.9,1.9)
9.75 (7.5,8.3)
2.71 (1.9,1.9)
+
11
MMSEnp
2.31 (−0.9,−0.9)
1.32 (1.6,1.7)
9.65 (7.5,8.2)
2.67 (1.9,1.9)
+
13
MAP
2.05 (0.6,0.6)
1.11 (1.5,1.5)
8.38 (6.5,6.3)
2.33 (1.7,1.7)
+
24
MMSE
1.79 (0.0,0.0)
0.96 (1.4,1.4)
7.93 (6.1,5.8)
2.03 (1.4,1.4)
+
33
ML
3.79 (−0.0,−0.0)
1.99 (3.0,3.0)
13.43 (13.2,11.7)
4.28 (3.0,3.0)
—
MAPnp
3.29 (−0.4,−0.4)
1.83 (2.7,2.6)
12.77 (12.6,11.5)
3.76 (2.7,2.7)
+
12
B
MMSEnp
3.10 (−1.2,−1.2)
1.68 (2.2,2.2)
10.99 (10.3,10.9)
3.52 (2.5,2.5)
+
18
MAP
2.56 (0.7,0.7)
1.34 (1.9,1.9)
8.95 (8.8,8.3)
2.89 (2.0,2.1)
+
33
MMSE
2.28 (0.0,0.0)
1.20 (1.8,1.8)
8.84 (8.6,8.1)
2.58 (1.8,1.8)
+
40
ML
1.99 (−0.0,−0.0)
1.13 (1.6,1.6)
8.40 (7.0,7.2)
2.29 (1.6,1.6)
—
MAPnp
1.81 (−0.2,−0.2)
1.05 (1.5,1.5)
7.00 (5.8,6.2)
2.09 (1.5,1.5)
+
9
C
MMSEnp
1.76 (−0.6,−0.6)
1.02 (1.3,1.3)
7.00 (5.8,6.2)
2.03 (1.4,1.5)
+
11
MAP
1.67 (0.4,0.4)
0.94 (1.3,1.3)
7.18 (5.5,5.7)
1.92 (1.4,1.4)
+
16
MMSE
1.48 (0.0,0.0)
0.82 (1.2,1.2)
6.95 (5.0,5.6)
1.70 (1.2,1.2)
+
26
ML
2.68 (0.0,0.0)
1.48 (2.2,2.2)
9.47 (7.8,9.0)
3.06 (2.2,2.2)
—
MAPnp
2.23 (−0.4,−0.4)
1.32 (1.8,1.8)
9.47 (7.4,7.9)
2.59 (1.8,1.8)
+
15
D
MMSEnp
2.35 (−1.2,−1.2)
1.32 (1.5,1.5)
9.45 (7.5,7.9)
2.69 (1.9,1.9)
+
12
MAP
1.38 (0.4,0.4)
0.76 (1.1,1.1)
5.78 (5.0,4.4)
1.57 (1.1,1.1)
+
49
MMSE
1.29 (0.0,0.0)
0.67 (1.0,1.0)
5.69 (4.6,4.2)
1.45 (1.0,1.0)
+
53
2
×
2
1.88 (0.0,0.0)
0.99 (1.5,1.5)
6.48 (6.4,6.2)
2.13 (1.5,1.5)
-
40
ML
1.33 (0.0,−0.0)
0.74 (1.1,1.1)
4.93 (3.8,4.2)
1.52 (1.1,1.1)
—
E
MAPnp
1.24 (−0.1,−0.1)
0.70 (1.0,1.0)
4.93 (3.8,4.2)
1.42 (1.0,1.0)
+
7
MMSEnp
1.20 (−0.3,−0.3)
0.68 (0.9,0.9)
4.93 (3.8,4.2)
1.38 (1.0,1.0)
+
9
MAP
1.22 (0.2,0.2)
0.67 (1.0,1.0)
4.74 (3.5,4.6)
1.39 (1.0,1.0)
+
9
MMSE
1.11 (0.0,0.0)
0.61 (0.9,0.9)
4.72 (3.4,4.6)
1.27 (0.9,0.9)
+
16Figure 8
Histogram of the residual radial errors of 2D parallel beam geometry for case (A) after estimation with the2×2 approach (radial RMSE: 4.24), ML (3.05),  MAPnp (2.71),  MMSEnp (2.67), MAP (2.33), and MMSE (2.03).The RMSEs of ML,MAPnp, MMSEnp, MAP, and MMSE for broadly varying parameters, such as the number of projections N (a), the standard deviation of the prior knowledge σX1=σX2 (b), and the standard deviation of the observation uncertainty σΔU (c). The constant parameters of each subfigure are set to the standard values of case (A). In the central subfigure MAPnp and MMSEnp are not presented since they are independent of σX1/2.
(a)
(b)
(c)Posterior distributions for the samex and full prior. (a) Case (B), using 2 projections. (b) Case (C), using 10 projections. (c) Case (D), using strong prior knowledge. (d) Case (E), using higher measurement accuracy.
(a)
(b)
(c)
(d)Based on  Table2  and  Figure 9  several characteristic observations can be made: at first, in all cases the estimation accuracy, which can be summarized with the RMSE, is improving mainly along the chain 2×2→ML→MAPnp→MMSEnp→MAP→MMSE.Second,  in order to compare the empirical bias of each estimator one has to focus on the mean values of the coordinate directionsx1 and x2 which are presented in brackets in the third column of Table 2. In theory, the 2×2 and the ML approaches are (asymptotically) unbiased, which was also observed in this empirical evaluation (i.e., a bias of 0.0 in each coordinate direction). MAP and MMSE estimators are typically biased estimators, which is certainly true in this application for the MAP estimation and partly true for the MMSE estimator. In detail, for the estimators with uniform priors, MAPnp and MMSEnp, a bias towards the center of the constraining circle can be recognized (e.g., negative bias values in each coordinate direction). This is plausible, since estimates outside the circle are prohibited with this prior and, on the other side, no information about the concentration of x at the border of the circle is introduced. When focusing on the full prior (constraining circle and the concentration at the border) with MAP, a bias in the opposite direction, towards the border, is introduced (e.g., see  Figure 7 for an illustration of this effect). Please note that the introduction of a significant bias in these three estimations, MAPnp, MMSEnp,  and MAP, still leads to better results compared to the unbiased 2×2 and ML estimations, which is a characteristic behavior in Bayesian estimation. Remarkably, the use of the Bayesian MMSE with the full prior leads to an empirically unbiased estimator in this application. This can be understood by realizing that the bias is in balance between the concentration of x at the border (leading to a bias towards the border) and the constraints of the border (leading to a bias towards the center of the circle) when averaging over the whole posterior distribution (e.g., see  Figure 7 again for an illustration). This allows to conclude that for the spatial estimation at constraining interfaces a highly accurate and practically unbiased estimator can be constructed with MMSE estimation.Third, compared to the standard parameter set, utilizingN=2 projections reduces the observation information and increases the RMSE. Specifically, ML is converging exactly to the 2×2 approach, increasing its RMSE from 3.05 to 4.28. Compared to the uniform prior in MAPnp and MMSEnp, the full prior in MAP and MMSE leads to robuster estimates with much less increasing RMSEs.Fourth, by utilizingN=10 projections, the performances of all estimators are improved. Essentially, the relative use of the prior knowledge is decreased since it is weighted less by the MAP and MMSE estimators, which, for example, can be seen in  (32).Fifth, when the sample distribution is getting narrower and simultaneously the corresponding prior PDF is utilized, the RMSEs of MAP and MMSE decrease significantly. The relative benefit from MAP to MMSE is decreased in this case, since, for the narrower prior, the influence of the truncation of the posterior by the boundary is reduced. On the other side,MAPnp and MMSEnp lead to RMSEs much closer to ML, since they introduce a strong bias (not knowing the concentration of x at the border). In this context, the increase in the bias is much stronger for MMSEnp compared to MAPnp, leading overall to even worse RMSE values.Sixth, if the measurement accuracy is increased, all estimators benefit from it significantly. Particularly, the benefit of using (partial) prior knowledge (inMAPnp, MMSEnp, MAP, and MMSE) compared to not using it (in 2×2 and ML) is decreased but still not negligible.


### 4.2. Cone Beam in 3D

In  Figure11 the estimation performances for this standard parameter set of the ML of  (39), the MAP of  (45), and the MMSE of  (46) are presented. Further, in  Table 3, the results of the three estimators are presented for all five cases (A)–(E), similar to the results of 2D parallel beam geometry.Table 3
Statistical values about the estimation results of 3D cone beam geometry for cases (A)–(E) withM= 10,000 sample points x each are presented for the three 3D estimators. The radial estimation error is presented at first and the coordinate estimation errors of (x1,x2,x3) are shown in brackets. The same formulas described as in the caption of Table 2 results are utilized.
Case
Method
Mean
st. deviation
abs. maximum
Radial RMSE
[
%
]
ML
2.52 (−0.0,−0.0,−0.0)
1.25 (1.8,1.8,1.1)
8.89 (7.7,7.1,4.7)
2.81 (1.8,1.8,1.1)
—
A
MAP
2.20 (0.8,0.8,0.3)
1.03 (1.4,1.4,1.0)
7.86 (6.1,6.9,4.1)
2.43 (1.6,1.6,1.0)
+
14
MMSE
1.77 (0.0,0.0,0.0)
0.78 (1.2,1.2,0.9)
5.99 (4.8,5.2,4.1)
1.93 (1.2,1.2,0.9)
+
31
ML
3.36 (0.1,0.1,0.0)
1.50 (2.3,2.3,1.7)
10.55 (9.5,9.1,7.2)
3.68 (2.3,2.3,1.7)
—
B
MAP
2.79 (0.9,0.9,0.6)
1.21 (1.6,1.6,1.4)
8.72 (7.9,8.3,5.5)
3.04 (1.9,1.9,1.5)
+
17
MMSE
2.27 (0.0,0.0,0.1)
1.00 (1.5,1.5,1.3)
7.80 (7.1,7.4,5.2)
2.48 (1.5,1.5,1.3)
+
33
ML
1.89 (0.0,0.0,0.0)
0.98 (1.4,1.4,0.8)
7.57 (5.5,6.9,2.9)
2.13 (1.4,1.4,0.8)
—
C
MAP
1.74 (0.6,0.6,0.2)
0.88 (1.1,1.1,0.7)
7.04 (4.8,5.5,3.0)
1.95 (1.3,1.3,0.8)
+
8
MMSE
1.42 (0.0,0.0,0.0)
0.68 (1.0,1.0,0.7)
5.18 (4.3,3.9,3.0)
1.57 (1.0,1.0,0.7)
+
36
ML
2.56 (0.0,0.0,0.0)
1.28 (1.9,1.9,1.1)
10.99 (7.4,8.7,4.7)
2.86 (1.9,1.9,1.1)
—
D
MAP
1.86 (0.9,0.9,0.5)
0.77 (0.9,0.9,0.8)
5.73 (5.7,5.3,3.6)
2.02 (1.3,1.3,0.9)
+
29
MMSE
1.31 (0.1,0.1,0.0)
0.59 (0.9,0.9,0.8)
5.00 (4.7,4.0,3.2)
1.43 (0.9,0.9,0.8)
+
50
ML
1.27 (−0.0,−0.0,−0.0)
0.63 (0.9,0.9,0.5)
4.78 (3.9,3.8,2.2)
1.41 (0.9,0.9,0.5)
—
E
MAP
1.21 (0.3,0.3,0.1)
0.60 (0.8,0.8,0.5)
4.55 (3.8,3.5,2.3)
1.35 (0.9,0.9,0.5)
+
4
MMSE
1.07 (0.0,0.0,0.0)
0.51 (0.8,0.8,0.5)
3.93 (3.0,2.8,2.0)
1.19 (0.8,0.8,0.5)
+
16Figure 11
Histogram of the residual radial errors of 3D cone beam geometry for case (A) after estimation with ML (radial RMSE: 2.81), MAP (2.43), and MMSE (1.93).Several observations can be made: at first, in general the same behavior as for the 2D parallel beam geometry can be made regarding the presented five cases, and, therefore, the mentioned statements remain valid; especially, the estimation improves in the order ML→ MAP → MMSE.Second, when focusing on the RMSE column, for all three methods the estimation of the third coordinatex3 is in general more accurate than the estimation of x1 and x2. This can be understood by the fact that x3 is almost directly observed in the cone beam geometry of  (34) and, therefore, is estimated with higher confidence.Third, the estimation ofx1 and x2 is slightly better in the 3D cone beam geometry than in the 2D parallel beam geometry. This is due to the fact that also in the third coordinate there is a slight information about the first two coordinates due to the scaling in cone beam geometry. Therefore, the information content about the first two coordinates is higher in 3D cone beam geometry, having two measurements per projection u for three unknowns x.Fourth, the calculation times on the test computer (Intel(R) Core2 Quad CPU, 4 Cores, 3.00 GHz, 4.00 GB RAM) for a single data point were the same for ML and MAP with in average1.9·10-3seconds utilizing the MATLAB implemented Nelder-Mead simplex method and for the MMSE implementation in average 5.5seconds utilizing 1/3 Simpson's rule for integration.  The optimization in ML and MAP was terminated if the difference of the estimated position between iterations was less than 10-6, and for the MMSE integration a 100×100×100 grid with a constant spatial step size of 0.2 was utilized. It could be observed that a further refinement of these criteria did not lead to practically better results. Several ways for a further speedup of the MMSE estimator can be implemented, for example, taking symmetries into account when utilizing Gaussians, more adaptive grid generation in the numerical integration, more effectively utilizing parallel programming and a full implementation in C compared to MATLAB.


## 4.1. Parallel Beam in 2D

For the standard parameters of case (A), in  Figure7 an example of the posterior PDF in 2D for a single sample point x is presented. In this plot on the right side a profile of the likelihood and the posterior function is presented, which demonstrates the general behavior: the width of the posterior is reduced compared to the likelihood function by the cost of introducing a shift (bias) of the maximum position by prior knowledge. Further, in the posterior profile the circular boundary is presented by the truncation of the profile. By incorporating the truncation into the prior distribution, the MMSE (which for simplicity one can think of as the expectation value of the posterior profile) is compensating for the bias of the maximum position in MAP.(a) Posterior density function for a specific estimation of case (A) and full prior. (b) Profile of the likelihood and posterior function atx=r·(cosπ/6,sinπ/6)T with r∈[0,12].
(a)
(b)In  Figure8 the estimation performance, for case (A) of the 2×2 approach of  (18), the ML of  (24), the MAP of  (32) and  (31), and the MMSE of  (33) (with noninformative prior and truncated Gaussian prior) are presented. The complete estimation results for all five cases (A)–(E) are presented in  Table 2.  In addition to  Table 2 the RMSEs of a broader variation of the essential parameters are illustrated in  Figure 9 by starting with case (A) and varying only the specified parameter of each subfigure. Examples for the posterior distributions  cases (A)–(E) (for the same specific sample x) are presented in  Figure 10, in order to see the qualitative differences due to the according changes in the parameter set.Table 2
Statistical values about the estimation results of 2D parallel beam geometry for cases (A)–(E) withM= 10,000 sample points x each are presented for the six 2D estimators. The radial estimation error is presented at first and the coordinate estimation errors of (x1,x2) are shown in brackets. In the following the formulas for the statistical values are presented for the radial error: the mean value=∥δ(u)∥-∥x∥¯=1/M∑j=1M‍∥δ(uj)∥-∥xj∥, the standard deviation=(1/(M-1)∑j=1M‍(∥δ(uj)∥-∥δ(u)∥¯-∥xj∥+∥x∥¯)2)1/2, the absolute maximum=maxj=1⋯M|∥δ(uj)∥-∥xj∥|, the Root Mean Square Error (RMSE) =(1/M∑j=1M‍(∥δ(uj)∥-∥xj∥2)1/2 and thegainof methodX relative to ML =(RMSEML-RMSEX)/RMSEML.
Case
Method
Mean
st. deviation
abs. maximum
Radial RMSE
[
%
]
2
×
2
3.77 (−0.0,0.0)
1.94 (3.0,3.0)
11.90 (11.1,10.9)
4.24 (3.0,3.0)
-
39
ML
2.67 (−0.0,0.0)
1.48 (2.2,2.2)
9.75 (9.2,8.3)
3.05 (2.2,2.2)
—
A
MAPnp
2.35 (−0.3,−0.3)
1.35 (1.9,1.9)
9.75 (7.5,8.3)
2.71 (1.9,1.9)
+
11
MMSEnp
2.31 (−0.9,−0.9)
1.32 (1.6,1.7)
9.65 (7.5,8.2)
2.67 (1.9,1.9)
+
13
MAP
2.05 (0.6,0.6)
1.11 (1.5,1.5)
8.38 (6.5,6.3)
2.33 (1.7,1.7)
+
24
MMSE
1.79 (0.0,0.0)
0.96 (1.4,1.4)
7.93 (6.1,5.8)
2.03 (1.4,1.4)
+
33
ML
3.79 (−0.0,−0.0)
1.99 (3.0,3.0)
13.43 (13.2,11.7)
4.28 (3.0,3.0)
—
MAPnp
3.29 (−0.4,−0.4)
1.83 (2.7,2.6)
12.77 (12.6,11.5)
3.76 (2.7,2.7)
+
12
B
MMSEnp
3.10 (−1.2,−1.2)
1.68 (2.2,2.2)
10.99 (10.3,10.9)
3.52 (2.5,2.5)
+
18
MAP
2.56 (0.7,0.7)
1.34 (1.9,1.9)
8.95 (8.8,8.3)
2.89 (2.0,2.1)
+
33
MMSE
2.28 (0.0,0.0)
1.20 (1.8,1.8)
8.84 (8.6,8.1)
2.58 (1.8,1.8)
+
40
ML
1.99 (−0.0,−0.0)
1.13 (1.6,1.6)
8.40 (7.0,7.2)
2.29 (1.6,1.6)
—
MAPnp
1.81 (−0.2,−0.2)
1.05 (1.5,1.5)
7.00 (5.8,6.2)
2.09 (1.5,1.5)
+
9
C
MMSEnp
1.76 (−0.6,−0.6)
1.02 (1.3,1.3)
7.00 (5.8,6.2)
2.03 (1.4,1.5)
+
11
MAP
1.67 (0.4,0.4)
0.94 (1.3,1.3)
7.18 (5.5,5.7)
1.92 (1.4,1.4)
+
16
MMSE
1.48 (0.0,0.0)
0.82 (1.2,1.2)
6.95 (5.0,5.6)
1.70 (1.2,1.2)
+
26
ML
2.68 (0.0,0.0)
1.48 (2.2,2.2)
9.47 (7.8,9.0)
3.06 (2.2,2.2)
—
MAPnp
2.23 (−0.4,−0.4)
1.32 (1.8,1.8)
9.47 (7.4,7.9)
2.59 (1.8,1.8)
+
15
D
MMSEnp
2.35 (−1.2,−1.2)
1.32 (1.5,1.5)
9.45 (7.5,7.9)
2.69 (1.9,1.9)
+
12
MAP
1.38 (0.4,0.4)
0.76 (1.1,1.1)
5.78 (5.0,4.4)
1.57 (1.1,1.1)
+
49
MMSE
1.29 (0.0,0.0)
0.67 (1.0,1.0)
5.69 (4.6,4.2)
1.45 (1.0,1.0)
+
53
2
×
2
1.88 (0.0,0.0)
0.99 (1.5,1.5)
6.48 (6.4,6.2)
2.13 (1.5,1.5)
-
40
ML
1.33 (0.0,−0.0)
0.74 (1.1,1.1)
4.93 (3.8,4.2)
1.52 (1.1,1.1)
—
E
MAPnp
1.24 (−0.1,−0.1)
0.70 (1.0,1.0)
4.93 (3.8,4.2)
1.42 (1.0,1.0)
+
7
MMSEnp
1.20 (−0.3,−0.3)
0.68 (0.9,0.9)
4.93 (3.8,4.2)
1.38 (1.0,1.0)
+
9
MAP
1.22 (0.2,0.2)
0.67 (1.0,1.0)
4.74 (3.5,4.6)
1.39 (1.0,1.0)
+
9
MMSE
1.11 (0.0,0.0)
0.61 (0.9,0.9)
4.72 (3.4,4.6)
1.27 (0.9,0.9)
+
16Figure 8
Histogram of the residual radial errors of 2D parallel beam geometry for case (A) after estimation with the2×2 approach (radial RMSE: 4.24), ML (3.05),  MAPnp (2.71),  MMSEnp (2.67), MAP (2.33), and MMSE (2.03).The RMSEs of ML,MAPnp, MMSEnp, MAP, and MMSE for broadly varying parameters, such as the number of projections N (a), the standard deviation of the prior knowledge σX1=σX2 (b), and the standard deviation of the observation uncertainty σΔU (c). The constant parameters of each subfigure are set to the standard values of case (A). In the central subfigure MAPnp and MMSEnp are not presented since they are independent of σX1/2.
(a)
(b)
(c)Posterior distributions for the samex and full prior. (a) Case (B), using 2 projections. (b) Case (C), using 10 projections. (c) Case (D), using strong prior knowledge. (d) Case (E), using higher measurement accuracy.
(a)
(b)
(c)
(d)Based on  Table2  and  Figure 9  several characteristic observations can be made: at first, in all cases the estimation accuracy, which can be summarized with the RMSE, is improving mainly along the chain 2×2→ML→MAPnp→MMSEnp→MAP→MMSE.Second,  in order to compare the empirical bias of each estimator one has to focus on the mean values of the coordinate directionsx1 and x2 which are presented in brackets in the third column of Table 2. In theory, the 2×2 and the ML approaches are (asymptotically) unbiased, which was also observed in this empirical evaluation (i.e., a bias of 0.0 in each coordinate direction). MAP and MMSE estimators are typically biased estimators, which is certainly true in this application for the MAP estimation and partly true for the MMSE estimator. In detail, for the estimators with uniform priors, MAPnp and MMSEnp, a bias towards the center of the constraining circle can be recognized (e.g., negative bias values in each coordinate direction). This is plausible, since estimates outside the circle are prohibited with this prior and, on the other side, no information about the concentration of x at the border of the circle is introduced. When focusing on the full prior (constraining circle and the concentration at the border) with MAP, a bias in the opposite direction, towards the border, is introduced (e.g., see  Figure 7 for an illustration of this effect). Please note that the introduction of a significant bias in these three estimations, MAPnp, MMSEnp,  and MAP, still leads to better results compared to the unbiased 2×2 and ML estimations, which is a characteristic behavior in Bayesian estimation. Remarkably, the use of the Bayesian MMSE with the full prior leads to an empirically unbiased estimator in this application. This can be understood by realizing that the bias is in balance between the concentration of x at the border (leading to a bias towards the border) and the constraints of the border (leading to a bias towards the center of the circle) when averaging over the whole posterior distribution (e.g., see  Figure 7 again for an illustration). This allows to conclude that for the spatial estimation at constraining interfaces a highly accurate and practically unbiased estimator can be constructed with MMSE estimation.Third, compared to the standard parameter set, utilizingN=2 projections reduces the observation information and increases the RMSE. Specifically, ML is converging exactly to the 2×2 approach, increasing its RMSE from 3.05 to 4.28. Compared to the uniform prior in MAPnp and MMSEnp, the full prior in MAP and MMSE leads to robuster estimates with much less increasing RMSEs.Fourth, by utilizingN=10 projections, the performances of all estimators are improved. Essentially, the relative use of the prior knowledge is decreased since it is weighted less by the MAP and MMSE estimators, which, for example, can be seen in  (32).Fifth, when the sample distribution is getting narrower and simultaneously the corresponding prior PDF is utilized, the RMSEs of MAP and MMSE decrease significantly. The relative benefit from MAP to MMSE is decreased in this case, since, for the narrower prior, the influence of the truncation of the posterior by the boundary is reduced. On the other side,MAPnp and MMSEnp lead to RMSEs much closer to ML, since they introduce a strong bias (not knowing the concentration of x at the border). In this context, the increase in the bias is much stronger for MMSEnp compared to MAPnp, leading overall to even worse RMSE values.Sixth, if the measurement accuracy is increased, all estimators benefit from it significantly. Particularly, the benefit of using (partial) prior knowledge (inMAPnp, MMSEnp, MAP, and MMSE) compared to not using it (in 2×2 and ML) is decreased but still not negligible.


## 4.2. Cone Beam in 3D

In  Figure11 the estimation performances for this standard parameter set of the ML of  (39), the MAP of  (45), and the MMSE of  (46) are presented. Further, in  Table 3, the results of the three estimators are presented for all five cases (A)–(E), similar to the results of 2D parallel beam geometry.Table 3
Statistical values about the estimation results of 3D cone beam geometry for cases (A)–(E) withM= 10,000 sample points x each are presented for the three 3D estimators. The radial estimation error is presented at first and the coordinate estimation errors of (x1,x2,x3) are shown in brackets. The same formulas described as in the caption of Table 2 results are utilized.
Case
Method
Mean
st. deviation
abs. maximum
Radial RMSE
[
%
]
ML
2.52 (−0.0,−0.0,−0.0)
1.25 (1.8,1.8,1.1)
8.89 (7.7,7.1,4.7)
2.81 (1.8,1.8,1.1)
—
A
MAP
2.20 (0.8,0.8,0.3)
1.03 (1.4,1.4,1.0)
7.86 (6.1,6.9,4.1)
2.43 (1.6,1.6,1.0)
+
14
MMSE
1.77 (0.0,0.0,0.0)
0.78 (1.2,1.2,0.9)
5.99 (4.8,5.2,4.1)
1.93 (1.2,1.2,0.9)
+
31
ML
3.36 (0.1,0.1,0.0)
1.50 (2.3,2.3,1.7)
10.55 (9.5,9.1,7.2)
3.68 (2.3,2.3,1.7)
—
B
MAP
2.79 (0.9,0.9,0.6)
1.21 (1.6,1.6,1.4)
8.72 (7.9,8.3,5.5)
3.04 (1.9,1.9,1.5)
+
17
MMSE
2.27 (0.0,0.0,0.1)
1.00 (1.5,1.5,1.3)
7.80 (7.1,7.4,5.2)
2.48 (1.5,1.5,1.3)
+
33
ML
1.89 (0.0,0.0,0.0)
0.98 (1.4,1.4,0.8)
7.57 (5.5,6.9,2.9)
2.13 (1.4,1.4,0.8)
—
C
MAP
1.74 (0.6,0.6,0.2)
0.88 (1.1,1.1,0.7)
7.04 (4.8,5.5,3.0)
1.95 (1.3,1.3,0.8)
+
8
MMSE
1.42 (0.0,0.0,0.0)
0.68 (1.0,1.0,0.7)
5.18 (4.3,3.9,3.0)
1.57 (1.0,1.0,0.7)
+
36
ML
2.56 (0.0,0.0,0.0)
1.28 (1.9,1.9,1.1)
10.99 (7.4,8.7,4.7)
2.86 (1.9,1.9,1.1)
—
D
MAP
1.86 (0.9,0.9,0.5)
0.77 (0.9,0.9,0.8)
5.73 (5.7,5.3,3.6)
2.02 (1.3,1.3,0.9)
+
29
MMSE
1.31 (0.1,0.1,0.0)
0.59 (0.9,0.9,0.8)
5.00 (4.7,4.0,3.2)
1.43 (0.9,0.9,0.8)
+
50
ML
1.27 (−0.0,−0.0,−0.0)
0.63 (0.9,0.9,0.5)
4.78 (3.9,3.8,2.2)
1.41 (0.9,0.9,0.5)
—
E
MAP
1.21 (0.3,0.3,0.1)
0.60 (0.8,0.8,0.5)
4.55 (3.8,3.5,2.3)
1.35 (0.9,0.9,0.5)
+
4
MMSE
1.07 (0.0,0.0,0.0)
0.51 (0.8,0.8,0.5)
3.93 (3.0,2.8,2.0)
1.19 (0.8,0.8,0.5)
+
16Figure 11
Histogram of the residual radial errors of 3D cone beam geometry for case (A) after estimation with ML (radial RMSE: 2.81), MAP (2.43), and MMSE (1.93).Several observations can be made: at first, in general the same behavior as for the 2D parallel beam geometry can be made regarding the presented five cases, and, therefore, the mentioned statements remain valid; especially, the estimation improves in the order ML→ MAP → MMSE.Second, when focusing on the RMSE column, for all three methods the estimation of the third coordinatex3 is in general more accurate than the estimation of x1 and x2. This can be understood by the fact that x3 is almost directly observed in the cone beam geometry of  (34) and, therefore, is estimated with higher confidence.Third, the estimation ofx1 and x2 is slightly better in the 3D cone beam geometry than in the 2D parallel beam geometry. This is due to the fact that also in the third coordinate there is a slight information about the first two coordinates due to the scaling in cone beam geometry. Therefore, the information content about the first two coordinates is higher in 3D cone beam geometry, having two measurements per projection u for three unknowns x.Fourth, the calculation times on the test computer (Intel(R) Core2 Quad CPU, 4 Cores, 3.00 GHz, 4.00 GB RAM) for a single data point were the same for ML and MAP with in average1.9·10-3seconds utilizing the MATLAB implemented Nelder-Mead simplex method and for the MMSE implementation in average 5.5seconds utilizing 1/3 Simpson's rule for integration.  The optimization in ML and MAP was terminated if the difference of the estimated position between iterations was less than 10-6, and for the MMSE integration a 100×100×100 grid with a constant spatial step size of 0.2 was utilized. It could be observed that a further refinement of these criteria did not lead to practically better results. Several ways for a further speedup of the MMSE estimator can be implemented, for example, taking symmetries into account when utilizing Gaussians, more adaptive grid generation in the numerical integration, more effectively utilizing parallel programming and a full implementation in C compared to MATLAB.


## 5. Discussion

We presented briefly the general Bayesian framework for point estimation and its application to the stochastic localization/triangulation problem. This compact presentation of the estimation methods is suitable for a broad group of applications. In the general introduction, we presented the foundations of Maximum Likelihood, Maximum A Posteriori, and Bayesian estimation applied to the Minimum Mean Square Error and how they are related to each other and other popular estimation methods.Several points can be made considering the specific application to stochastic localization/triangulation, which hold true in general for Bayesian approaches as well: as seen in the proposed application, there is a tradeoff between benefit and effort. The more the  information is contained in observations (= quality of observation, e.g., more projections or higher quality measurement/detection), the smaller the benefit is due to prior knowledge. The other way round is true as well: the less the information in observations, the more important is the prior knowledge. With the Bayesian estimation (MAPnp, MMSEnp, MAP, or MMSE), these two probabilistic sources of information can be naturally included in the model.The introduction of aprior can be utilized to benefit the estimation in two ways: it can place more weight on possible values of locations x and effectively restrict ΩX by setting the prior probabilities to zeros outside of the feasibility region (boundaries). If the constrained region of interest for x is a convex region (as in the example), then the MMSE and the MAP will derive estimates in the constrained region, and, therefore, this imposes an effective constraint. However, if the region is nonconvex, only the MAP will still lead to an estimate within the restricted region, while MMSE will lead to an estimate in the convex hull of the constrained region (possibly outside the constrained region). In a broader context, the prior can also be interpreted as a regularizing influence to the presented inverse problem, such as the Tikhonov regularization as in MAP or in a more interweaved form in MMSE.It is pointed out that in this simulation study the standard deviations of the probabilistic knowledge about the prior knowledge (σXi for i={1,2,3}) and the observation uncertainties (σΔU) are known before estimation. In some practical applications these parameters are not known and have to be estimated or assumed, which would introduce additional uncertainties to the presented estimation methods. The investigation of these additional uncertainties is not part of this study, and it is recommended to extend the presented methods for this purpose utilizing hierarchical Bayesian models [13].Focusing on the presented 2D study, the following statements can be formulated: if no prior knowledge is available, utilizing ML with 5 or 10 projections compared to a single2×2 with 2 projections decreases effectively the measurement uncertainties and increases the estimation accuracy by about 40%. If only information about the constraints of x is available but no information about its possible distribution (e.g., truncated Gaussian), MAPnp should be preferred to MMSEnp since it achieves a similar accuracy (case (A): 11% increase for MAPnp relative to ML compared to 13% for MMSEnp) with less bias and much faster calculation times. In contrast to this, if both information are available, information about the constraints and concentration of x, then utilizing MMSE should be preferred to MAP since it allows finding a better compromise close to the surfaces between constraints and concentration. This is due to the fact that MAP neglects the constraints totally as long as the estimates are inside the circle, leading to worse results in average (case (A): 33% increase for MMSE relative to ML compared to 24% for MAP).  In general, it is shown in this application that with the utilization of the MMSE for spatial coordinate estimation at constraining interfaces it is possible to construct an empirically unbiased estimator with high accuracy.Further testing of these methods in medical imaging or other practical applications of stochastic localization at boundaries is highly encouraged. This paper summarizes the basic theory for this purpose and provides clearly interpretable results in the simulation studies. However, in each practical application the true benefit of the estimation methods has to be evaluated separately or individually. For instance, in radiotherapy due to differences in patient anatomy, the problem of breast tumor localization is dramatically different from prostate localization. For breast the boundary is the breast surface, and, for prostate, hard constraints can arise from the bony anatomy or an endorectal balloon which fixes the anatomy differently.The presented estimation with constraints incorporated into the prior is very general and not limited to linear acquisition geometries, as demonstrated for the nonlinear 3D case. Further, the results of the 3D study show that the scaling of the 3D cone beam geometry increases the estimation accuracy for the first two coordinates ofx compared to the parallel beam geometry in 2D. Eventually, ML and MAP estimations lead to fast optimization problems and the MMSE estimator to a more challenging high-dimensional integration. In practical situations it depends on the specific application, such as real-time imaging, to find a suitable compromise between computational effort and accuracy.


## 6. Conclusion

With the 2D parallel beam and 3D cone beam geometry examples it is illustrated that the general framework of Bayesian estimation can be applied to a variety of estimation problems even when no explicit formulas can be derived, such as for nonlinear/non-Gaussian distributions of the uncertainties or prior knowledge. In this context, it is demonstrated that constraints about the solution space can be efficiently incorporated into the prior for stochastic localization by utilizing Maximum A Posteriori and Minimum Mean Square Error estimation.



---

*Source: 960421-2013-05-07.xml*