# Optimal estimation

In applied statistics, optimal estimation is a regularized matrix inverse method based on Bayes theorem. It is used very commonly in the geosciences, particularly for atmospheric sounding. A matrix inverse problem looks like this:

$\mathbf {A} {\vec {x}}={\vec {y}}$ The essential concept is to transform the matrix, A, into a conditional probability and the variables, ${\vec {x}}$ and ${\vec {y}}$ into probability distributions by assuming Gaussian statistics and using empirically-determined covariance matrices.

## Derivation

Typically, one expects the statistics of most measurements to be Gaussian. So for example for $P({\vec {y}}|{\vec {x}})$ , we can write:

$P({\vec {y}}|{\vec {x}})={\frac {1}{(2\pi )^{mn/2}|{\boldsymbol {S_{y}}}|}}\exp \left[-{\frac {1}{2}}({\boldsymbol {A}}{\vec {x}}-{\vec {y}})^{T}{\boldsymbol {S_{y}}}^{-1}({\boldsymbol {A}}{\vec {x}}-{\vec {y}})\right]$ $P({\vec {x}})={\frac {1}{(2\pi )^{m/2}|{\boldsymbol {S_{x_{a}}}}|}}\exp \left[-{\frac {1}{2}}({\vec {x}}-{\widehat {x_{a}}})^{T}{\boldsymbol {S_{x_{a}}}}^{-1}({\vec {x}}-{\widehat {x_{a}}})\right]$ The nice thing about the Gaussian distributions is that only two parameters are needed to describe them and so the whole problem can be converted once again to matrices. Assuming that $P({\vec {x}}|{\vec {y}})$ takes the following form:

$P({\vec {x}}|{\vec {y}})={\frac {1}{(2\pi )^{mn/2}|{\boldsymbol {S_{x}}}|}}\exp \left[-{\frac {1}{2}}({\vec {x}}-{\widehat {x}})^{T}{\boldsymbol {S_{x}}}^{-1}({\vec {x}}-{\widehat {x}})\right]$ ${\boldsymbol {S_{x}}}=({\boldsymbol {A}}^{T}{\boldsymbol {S_{y}^{-1}}}{\boldsymbol {A}}+{\boldsymbol {S_{x_{a}}^{-1}}})^{-1}$ ${\widehat {x}}={\widehat {x_{a}}}+{\boldsymbol {S_{x}}}{\boldsymbol {A}}^{T}{\boldsymbol {S_{y}}}^{-1}({\vec {y}}-{\boldsymbol {A}}{\widehat {x_{a}}})$ Because we are using Gaussians, the expected value is equivalent to the maximum likely value, and so this is also a form of maximum likelihood estimation.

Typically with optimal estimation, in addition to the vector of retrieved quantities, one extra matrix is returned along with the covariance matrix. This is sometimes called the resolution matrix or the averaging kernel and is calculated as follows:

${\boldsymbol {R}}=({\boldsymbol {A}}^{T}{\boldsymbol {S_{y}}}^{-1}{\boldsymbol {A}}+{\boldsymbol {S_{x_{a}}}}^{-1})^{-1}{\boldsymbol {A}}^{T}{\boldsymbol {S_{y}}}^{-1}{\boldsymbol {A}}$ This tells us, for a given element of the retrieved vector, how much of the other elements of the vector are mixed in. In the case of a retrieval of profile information, it typical indicates the altitude resolution for a given altitude. For instance if the resolution vectors for all the altitudes contain non-zero elements (to a numerical tolerance) in their four nearest neighbours, then the altitude resolution is only one fourth that of the actual grid size.