Point estimation

The method of maximum likelihood, due to R.A. Fisher, is arguably the most important general method for estimating the parameters ${\displaystyle \theta _{1},\dots ,\theta _{p}}$ of a parametric model.^[12] According to this method, the best estimate for the unknown model parameter is the one that maximizes the probability (or probability density) of observing the data that has been observed. This probability as a function of the model parameters is called the likelihood, giving the method its name.^[13]

In mathematical terms, the methods works as follows: Let ${\displaystyle X=(X_{1},\dots ,X_{n})}$ denote a random data sample with joint probability density function or probability mass function ${\displaystyle f_{\theta }}$ depending on the vector ${\displaystyle \theta =(\theta _{1},\dots ,\theta _{p})}$ of model parameters. The function ${\displaystyle \theta \mapsto f_{\theta }(X)}$ is called the likelihood function, often denoted by ${\displaystyle L}$. Assuming that true model parameter lies in a set ${\displaystyle \Theta \subset \mathbb {R} ^{p}}$, a maximum likelihood estimator ${\displaystyle {\hat {\theta }}_{\mathrm {MLE} }}$ fulfills

$${\displaystyle L({\hat {\theta }}_{\mathrm {MLE} })=\max _{\theta \in \Theta }L(\theta )=\max _{\theta \in \Theta }f_{\theta }(X).}$$

In practice, the likelihood function is often differentiable, in which case ${\displaystyle {\hat {\theta }}_{\mathrm {MLE} }}$ is a solution of the equations^[14]

$${\displaystyle {\frac {\partial L(\theta )}{\partial \theta _{i}}}=0,\quad i=1,\dots ,k.}$$

Under certain assumptions on the likelihood, the MLE is strongly consistent, asymptotically efficient and asymptotically normal.

The method of moments is one of the oldest methods of estimation, which was first introduced by Karl Pearson in 1895.^[15][16][17] This method uses the functional relationship between quantities of interest of a distribution (e.g. the variance or model parameters) and the moments of that distribution, leading to a set of equations that are solved for the aforementioned quantities. Motivated by the law of large numbers, the moments of the distribution are estimated by the sample moments, that is, the arithmetic averages of powers of the data values.^[18] The method stands out for its simplicity and wide domain of applicatibility, but often do not yield the best estimators and sometimes even no estimator at all.^[19] There also exists generalizations of the method of moments.

Mathematically speaking, for a distribution ${\displaystyle \mathbb {P} _{\theta }}$ depending on some vector of parameters ${\displaystyle \theta =(\theta _{1},\dots ,\theta _{k})}$, the equations of the first ${\displaystyle k}$ moments of ${\displaystyle X\sim \mathbb {P} _{\theta }}$ are established:

$${\displaystyle {\begin{aligned}&m_{1}=\mathbb {E} [X]=f_{1}(\theta _{1},\dots ,\theta _{k}),\\&\dots \\&m_{k}=\mathbb {E} [X^{k}]=f_{k}(\theta _{1},\dots ,\theta _{k}).\end{aligned}}}$$

If these equations have a solution, then there exists functions ${\displaystyle g_{1},\dots ,g_{k}}$ such that the parameters can be expressed by the moments via

$${\displaystyle {\begin{aligned}&\theta _{1}=g_{1}(m_{1},\dots ,m_{k}),\\&\dots \\&\theta _{k}=g_{k}(m_{1},\dots ,m_{k}).\end{aligned}}}$$

Given a random sample ${\displaystyle X_{1},\dots ,X_{n}}$ from ${\displaystyle \mathbb {P} _{\theta }}$, the moments are then replaced by the sample moments ${\displaystyle {\hat {m}}_{j}={\frac {1}{n}}\sum _{i=1}^{n}X_{i}^{j}}$, ${\displaystyle j=1,\dots ,k}$, yielding the moment estimators:

$${\displaystyle {\begin{aligned}&{\hat {\theta }}_{1}=g_{1}({\hat {m}}_{1},\dots ,{\hat {m}}_{k}),\\&\dots \\&{\hat {\theta }}_{k}=g_{k}({\hat {m}}_{1},\dots ,{\hat {m}}_{k}).\end{aligned}}}$$

Under relatively low assumptions on the solution functions ${\displaystyle g_{1},\dots ,g_{k}}$ and the moments of ${\displaystyle \mathbb {P} _{\theta }}$, the moment estimators are strongly consistent and asymptotically normal.

The Lehmann–Scheffé theorem, although not an estimation method in itself, motivates two methods for constructing UMVU estimators in models, where a complete, sufficient statistic ${\displaystyle T}$ is available.^[20][21] This applies, for example, to exponential families, where the natural sufficient statistic is always complete.

Let ${\displaystyle X=(X_{1},\dots ,X_{n})}$ be a vector of random samples from some distribution ${\displaystyle \mathbb {P} _{\theta }}$ depending on a parameter ${\displaystyle \theta }$.

Method 1: For the estimation of ${\displaystyle \theta }$, there exists at most one function ${\displaystyle g}$ that solves the equation $${\displaystyle \mathbb {E} _{\theta }[g(T(X))]=\theta ,}$$ in which case the estimator ${\displaystyle {\hat {\theta }}=g(T)}$ is an UMVU estimator for ${\displaystyle \theta }$.

Method 2: Given any unbiased estimator ${\displaystyle {\tilde {\theta }}}$ of ${\displaystyle \theta }$, the estimator $${\displaystyle {\hat {\theta }}(X)=\mathbb {E} _{\theta }[{\tilde {\theta }}(X)\mid T(X)],}$$ that is, the conditional expectation of ${\displaystyle {\tilde {\theta }}}$ given ${\displaystyle T}$, is the unique UMVU estimator. In some cases, computing the conditional expectation of a suitably chosen unbiased estimator given ${\displaystyle T}$ might be easier than solving the equation from method 1.

Score matching is a rather new method for parameter estimation that is less intuitive than MLE, but is more advantageous for complicated distributions with a large number of parameters.^[22]

The data is assumed to be drawn from a distribution with probability density ${\displaystyle f(\mathbf {x} ,{\boldsymbol {\theta }})}$, ${\displaystyle \mathbf {x} \in \mathbb {R} ^{d}}$, parameterized by ${\displaystyle {\boldsymbol {\theta }}}$. The key quantity, that is needed for the score matching estimation, is the score${\displaystyle \nabla _{\mathbf {x} }\log f(\mathbf {x} ,{\boldsymbol {\theta }})}$. The true parameter (vector) ${\displaystyle {\boldsymbol {\theta }}}$ is characterized by$${\displaystyle {\begin{aligned}{\boldsymbol {\theta }}&=\mathrm {arg} \min _{{\boldsymbol {\theta }}'}{\frac {1}{2}}\int _{\mathbb {R} ^{d}}\lVert \nabla _{\mathbf {y} }\log f(\mathbf {y} ,{\boldsymbol {\theta }}')-\nabla _{\mathbf {y} }\log f(\mathbf {y} ,{\boldsymbol {\theta }})\rVert ^{2}f(\mathbf {y} ,{\boldsymbol {\theta }})\mathrm {d} \mathbf {y} \\&=\mathrm {arg} \min _{{\boldsymbol {\theta }}'}\sum _{i=1}^{d}\int _{\mathbb {R} ^{d}}\left(\partial _{y_{i}}^{2}\log f(\mathbf {y} ,{\boldsymbol {\theta }}')+{\frac {1}{2}}(\partial _{y_{i}}\log f(\mathbf {y} ,{\boldsymbol {\theta }}')^{2}\right)f(\mathbf {y} ,{\boldsymbol {\theta }})\mathrm {d} \mathbf {y} .\end{aligned}}}$$In line with the second equality, for observations ${\displaystyle \mathbf {x} _{1},\dots ,\mathbf {x} _{n}}$, the score matching estimate ${\displaystyle {\hat {\boldsymbol {\theta }}}}$ is defined by$${\displaystyle {\hat {\boldsymbol {\theta }}}=\mathrm {arg} \min _{{\boldsymbol {\theta }}'}\sum _{j=1}^{n}\sum _{i=1}^{d}\left(\partial _{x_{i}}^{2}\log f(\mathbf {x} _{j},{\boldsymbol {\theta }}')+{\frac {1}{2}}(\partial _{x_{i}}\log f(\mathbf {x} _{j},{\boldsymbol {\theta }}')^{2}\right).}$$The main advantage of score matching estimation compared to the maximum likelihood estimation is the following. Typically, the probability density is of the form$${\displaystyle f(\mathbf {x} ,{\boldsymbol {\theta }})={\frac {1}{Z({\boldsymbol {\theta }})}}q(\mathbf {x} ,{\boldsymbol {\theta }}),}$$with an unnormalized density ${\displaystyle q(\mathbf {x} ,{\boldsymbol {\theta }})}$ and a normalization constant${\displaystyle Z({\boldsymbol {\theta }})}$. The computation of latter is often the bottleneck of maximum likelihood estimation, especially if the dependence on ${\displaystyle {\boldsymbol {\theta }}}$ is non-analytic and the parameter vector is high-dimensional. By taking the derivatives of the logarithm of ${\displaystyle f}$, ${\displaystyle Z({\boldsymbol {\theta }})}$ does not occur in the score matching estimation procedure, thus eliminating this issue.^[22]

Additionally, under mild restrictions on ${\displaystyle f(\mathbf {x} ,{\boldsymbol {\theta }})}$, the score matching estimate can be shown to be weakly consistent.

The method of least squares is a parameter estimation procedure for regression, where observations ${\displaystyle Y_{1},\dots ,Y_{n}}$ need to be expressed as functions of the covariates${\displaystyle X_{1},\dots ,X_{n}\in \mathbb {R} ^{p}}$ via

$${\displaystyle Y_{i}=f_{\beta }(X_{i})+\varepsilon _{i}}$$

for some function ${\displaystyle f_{\beta }}$ depending on a parameter vector ${\displaystyle \beta \in \mathbb {R} ^{p+1}}$ and a random variable ${\displaystyle \varepsilon _{i}}$ accounting for random noise in the data. The ordinary least squares estimator (OLS) ${\displaystyle {\hat {\beta }}_{\mathrm {OLS} }}$ assumes that ${\displaystyle f_{\beta }}$ is a linear function of ${\displaystyle \beta }$ and estimates the parameter by minimizing of the sum of squared differences, called residuals:

$${\displaystyle {\hat {\beta }}_{\mathrm {OLS} }\in \mathrm {arg} \min _{\beta \in \mathbb {R} ^{p}}\sum _{i=1}^{n}\left(Y_{i}-\beta _{0}-X_{i,1}\beta _{1}-\dots -X_{i,p}\beta _{p}\right)^{2}.}$$

In fact, if the noise variables ${\displaystyle \varepsilon _{i}}$ are independent and identically normally distributed with mean zero, then the OLS coincides with the maximum likelihood estimator.

If the rank of the matrix of covariates ${\displaystyle X=(X_{1},\dots ,X_{n})^{\top }\in \mathbb {R} ^{n\times p}}$ is equal to ${\displaystyle p}$, then the OLS is uniquely given by

$${\displaystyle {\hat {\beta }}_{\mathrm {OLS} }=(X^{\top }X)^{-1}X^{\top }Y,}$$

and it is the best linear unbiased estimator (BLUE) by the Gauß–Markov theorem. If, as above, the noise variables are independent and identically normally distributed, then the OLS is even UMVU for the entire class of unbiased estimators.^[23]

There are multiple variants of least squares estimation:

• weighted least squares (WLS) estimation: The residuals are weighted with different weighting factors ${\displaystyle w_{i}}$.
• generalized least squares (GLS) estimation: The noise variables are allowed to be correlated.
• nonlinear least squares estimation: The function ${\displaystyle f_{\beta }}$ can depend non-linearly on the parameter ${\displaystyle \beta }$.
• regularized least squares estimation: The parameter is chosen such that the sum of squared residuals plus a regularization term is minimized. Typical choices for the regularization are the euclidian norm ${\displaystyle \lVert \beta \rVert _{2}^{2}}$, leading to Ridge regression, or the 1-norm ${\displaystyle \lVert \beta \rVert _{1}}$, leading to Lasso regression.

Given a vector of random observations ${\displaystyle X=(X_{1},\dots ,X_{n})}$ and a loss function${\displaystyle L}$, a Bayes estimator is any estimator ${\displaystyle {\hat {\theta }}(X)}$ that minimizes the Bayes risk:^[24]$${\displaystyle R_{\pi }({\hat {\theta }})=\int _{\Theta }\mathbb {E} [L(\theta ,{\hat {\theta }}(X))]\pi (\theta )\mathrm {d} \theta =\mathbb {E} \left[\int _{\Theta }L{\big (}\theta ,{\hat {\theta }}(X){\big )}p(\theta \mid X)\mathrm {d} \theta \right].}$$

Obviously, the risk is minimized when the integral under the expectation is minimized as a function of the observations. Thus, a Bayes estimator is equivalently given as a minimizer of the posterior loss:

$${\displaystyle {\hat {\theta }}(X)\in \mathrm {arg} \min _{y\in \mathbb {R} ^{n}}\int _{\Theta }L{\big (}\theta ,y{\big )}p(\theta \mid X)\mathrm {d} \theta .}$$

The most common choices for ${\displaystyle L}$ are the following:

• ${\displaystyle L(\theta ,y)=(\theta -y)^{2}}$: In this case, the posterior loss is minimized by the posterior mean${\displaystyle {\hat {\theta }}(X)=\int _{\Theta }\theta \,p(\theta \mid X)\mathrm {d} \theta }$, as observed by Gauss.^[25]

• ${\displaystyle L(\theta ,y)=|\theta -y|}$: Here, the posterior loss is minimized by the median of the posterior distribution, as observed by Laplace.^[25]

By Wald's theorem, unique Bayesian estimators are admissible.^[26][27] Moreover, under certain assumptions on the prior and ${\displaystyle p_{\theta }}$, the posterior mean estimator can be shown to be consistent, asymptotically normal and asymptotically efficient.^[28]

The maximum a posteriori estimator ${\displaystyle {\hat {\theta }}_{\mathrm {MAP} }}$ is a maximizer of the posterior distribution:

$${\displaystyle {\hat {\theta }}_{\mathrm {MAP} }(X)\in \mathrm {arg} \max _{\theta \in \Theta }p(\theta \mid X).}$$

In a sense, the MAP estimator is a Bayesian analog to the maximum likelihood estimator. In fact, for a uniform prior distribution, the MAP estimator coincides with the MLE.

The MAP estimator has good asymptotic properties, even for many difficult problems, on which the maximum-likelihood estimator has difficulties. For regular problems, where the maximum-likelihood estimator is consistent, the maximum-likelihood estimator ultimately agrees with the MAP estimator.^[29][26][30]

The Minimum Message Length (MML) point estimator is based in Bayesian information theory and is not so directly related to the posterior distribution.

Special cases of Bayesian filters are important:

• Kalman filter
• Wiener filter

Several methods of computational statistics have close connections with Bayesian analysis:

• particle filter
• Markov chain Monte Carlo (MCMC)