Statistical inference is about drawing conclusions from data. This process often involves estimating some parameters of interest. In a parametric model, the parameters of interest are often the parameters of the corresponding parametric family. For a parametric model, there are three common tasks in statistical inference – estimating the underlying parameter, providing an interval inference about the underlying parameter, and testing if the underlying parameter satisfies certain conditions.
In statistics, there are two major paradigms for making inference – the Frequentist paradigm and the Bayesian paradigm1. We will talk about their principles of estimation, interval inference, and testing.
Note that there is no right or wrong about each paradigm – they are just different ways of making arguments. Each paradigm is a self-consistent way of making logical arguments and has its own advantages and limitations.
The Frequentist approach is the paradigm we learn from Statistics 101. It interprets probability as the long-term frequency. In the Frequentist approach, the parameter of interest is a fixed and unknown number.
In a parametric model, we often estimate the parameter of interest using the so-called maximum likelihood estimator (MLE). The idea is very simple. Suppose we observe only one observation \(X\) from a PDF/PMF \(p(x)\). The parametric model assumes that such a PDF/PMF can be written as \(p(x) = p(x;\theta)\), where \(\theta\) is the parameter of the model (\(\theta\) is often the parameter of interest) inside a parameter space \(\Theta\) (\(\theta\in \Theta\)). The idea of MLE is to ask the following question: given the observation \(X\), which \(\theta\) is the most likely parameter that generates \(X\)? To answer this question, we can vary \(\theta\) and examine the value of \(p(X;\theta)\).
Because we are treating \(X\) as fixed and \(\theta\) being something that we want to optimize, we can view the problem as finding the best \(\theta\) such that the likelihood function \(L(\theta|X) = p(X;\theta)\) is maximized. The MLE uses the \(\theta\) that maximizes the likelihood value. Namely, \[\hat{\theta}_{MLE} = {\sf argmax}_\theta L(\theta|X).\]
When we have multiple observations \(X_1,\cdots, X_n\), the likelihood function can be defined in a similar way – we use the joint PDF/PMF to define the likelihood function. Let \(p(x_1,\cdots,x_n;\theta)\) be the joint PDF/PMF. Then the likelihood function is \[L_n(\theta) = L(\theta|X_1,\cdots,X_n) = p(X_1,\cdots,X_n;\theta).\] Note that when we assume IID observations, \[L_n(\theta) = \prod_{i=1}^n L(\theta|X_i)=\prod_{i=1}^n p(X_i;\theta).\]
In many cases, instead of using the likelihood function, we often work with the log-likelihood function \[\ell_n(\theta) = \log L_n(\theta).\] Because taking the logarithm does not change the maximizer of a function, the maximizer of the log-likelihood function is the same as the maximizer of the likelihood function. There are both computational and mathematical advantages of using a log-likelihood function over a likelihood function. To see this, we consider the case of an IID sample. Computationally, the likelihood function often has a very small value due to the product form of PDF/PMFs. So it is very likely that the number is too small, making computation very challenging. Mathematically, when we take the log of the likelihood function, the product of PDF/PMFs becomes an additive form \[\ell_n(\theta) = \log L_n(\theta) = \sum_{i=1}^n \log p(X_i;\theta) .\] Under the IID assumption, each \(\log p(X_i;\theta)\) is an IID random variable so the central limit theorem and the Law of Large Numbers can be applied to the average, making it possible to analyze its asymptotic behavior.
Since under the IID assumption we have many advantages, we will assume IID from now on. Because MLE finds the maximum of \(\ell_n(\theta)\), a common trick to find the MLE is to study the gradient of the log-likelihood function, which is also known as the score function: \[S_n(\theta) = \frac{\partial}{\partial\theta} \ell_n(\theta) = \sum_{i=1}^n s(\theta|X_i),\] where \(s(\theta|X_i) = \frac{\partial}{\partial\theta} \ell(\theta|X_i)= \frac{\partial}{\partial\theta} \log p(X_i;\theta)\). Under suitable conditions, the MLE satisfies the score equation: \[S_n(\hat{\theta}_{MLE}) = 0.\] Note that if there is more than one parameter, say \(\theta\in \mathbb{R}^p\), the score equation will be a system of \(p\) equations.
Because the MLE is at the maximum point of the likelihood function, the curvature of the likelihood function around the maximum will determine its stability. To measure the curvature, we use the Fisher information matrix: \[I_n(\theta) = -\mathbb{E}\left[\frac{\partial^2}{\partial\theta \partial\theta^T} \ell_n(\theta)\right] = n\cdot I_1(\theta) = n\cdot -\mathbb{E}\left[\frac{\partial^2}{\partial\theta \partial\theta^T} p(X_1;\theta)\right].\]
If the data is generated from a PDF/PMF \(p(x;\theta_0)\) and some regularity conditions are satisfied, \[\begin{align*} \mathbb{E}(S_n(\theta_0)) &= 0,\\ I_1(\theta_0) &= \mathbb{E}(S_1(\theta_0)S^T_1(\theta_0)). \end{align*}\] Moreover, \[\sqrt{n}\left(\hat{\theta}_{MLE} - \theta_0\right)\overset{D}{\rightarrow} N(0, I^{-1}_1(\theta_0)).\] Namely, the MLE is asymptotically normally distributed around the true parameter \(\theta_0\) and the covariance is determined by the Fisher information matrix. Note that the asymptotic normality also implies that \(\hat{\theta}_{MLE}-\theta_0 \overset{P}{\rightarrow}0\).
Example 1: Binomial Distribution. Assume that we obtain a single observation \(Y\sim {\sf Bin}(n,p)\), and we assume that \(n\) is known. The goal is to estimate \(p\). The log-likelihood function is \[\ell(p) = Y\log p + (n-Y) \log (1-p) + C_n(Y),\] where \(C_n (Y)= \log { n\choose Y}\) is independent of \(p\). The score function is \[S(p) = \frac{Y}{p} -\frac{n-Y}{1-p}\] so solving the score equation gives us \(\hat{p}_{MLE} = \frac{Y}{n}\). Moreover, the Fisher information is \[I(p) =\mathbb{E}\left\{\frac{\partial}{\partial p} S(p)\right\} = -\frac{\mathbb{E}(Y)}{p^2} - \frac{n-\mathbb{E}(Y)}{(1-p)^2} = \frac{n}{p(1-p)}.\]
Example 2: Multinomial Distribution. Let \(X_1,\cdots,X_n\) be IID from a multinomial distribution such that \(P(X_1=j) = p_j\) for \(j=1,\cdots, s\) and \(\sum_{j=1}^s p_j = 1\). Note that the parameter space is \(\Theta = \{(p_1,\cdots, p_s): 0\leq p_j, \sum_{j=1}^s p_j=1\}\). By setting \(N_j = \sum_{i=1}^n I(X_i = j)\) for each \(j=1,\cdots, s\), we obtain the random vector \((N_1,\cdots, N_s)\sim {\sf Multinomial}(n,{p}),\) where \({p} = (p_1,\cdots, p_s)\). The parameters of interest are \(p_1,\cdots, p_s\). In this case, the likelihood function is \[L_n(p_1,\cdots, p_s) = \frac{n!}{N_1!\cdots N_s!} p_1^{N_1}\cdots p_s^{N_s}\] and the log-likelihood function is \[\ell_n(p_1,\cdots, p_s) = \sum_{j=1}^s N_j\log p_j + C_n,\] where \(C_n\) is independent of \(p\). Note that naively computing the score function and setting it to 0 will not grant us a solution (think about why) because we do not use the constraint of the parameter space – the parameters sum to 1. To use this constraint in our analysis, we consider adding Lagrange multipliers and optimizing it: \[F(p, \lambda) = \sum_{j=1}^s N_j\log p_j + \lambda\left(1-\sum_{j=1}^s p_j\right).\] Differentiating this function with respect to \(p_1,\cdots, p_s,\) and \(\lambda\) and setting them to 0 gives \[\begin{align*} \frac{\partial F}{\partial p_j} = \frac{N_j}{p_j} - \lambda 0 \Rightarrow N_j = \lambda \hat{p}_{MLE,j} \end{align*}\] and \(1-\sum_{j=1}^s p_j = 0\). Thus, \(n = \sum_{j=1}^s N_j = \lambda \sum_{j=1}^p = \lambda\) so \(\hat{p}_{MLE, j} =\frac{N_j}{n}\).
In some analysis, we not only want to have just a point estimate of the parameter of interest, but also want to use an interval to infer the parameter of interest. And we also want to assign a level to this interval to describe how ‘accurate’ this interval is. Note that here the concept of accuracy is not well-defined – we will talk about it later. Ideally, given an accuracy level, we want an interval as small as possible.
Frequentists and Bayesians define accuracy differently so their constructions of intervals are also different. In short, Frequentists define accuracy as the long-term frequency coverage of the underlying true parameter of interest whereas Bayesians define accuracy in terms of covering posterior probability. In this section, we will talk about the Frequentist approach and the interval is known as the confidence interval. The accuracy that Frequentists use is called the confidence level.
Formally, given a confidence level \(1-\alpha\), a confidence interval of \(\theta_0\) is a random interval \(C_{n,\alpha}\) that can be constructed solely from the data (i.e., can be constructed using \(X_1,\cdots, X_n\)) such that \[P(\theta\in C_{n,\alpha}) \geq 1-\alpha +o(1).\] Beware, what is random is not \(\theta\) but the interval \(C_{n,\alpha}\). The quantity \(P(\theta\in C_{n,\alpha})\) is also called the (Frequentist) coverage. Note that we allow the coverage to be asymptotically \(1-\alpha\); when there is no \(o(1)\) term, we will say that the confidence interval has finite sample coverage. A confidence interval with the above property is also called an (asymptotically) valid confidence interval.
Normal confidence interval. A traditional approach to constructing a confidence interval of \(\theta_0\) is based on the asymptotic normality of the MLE: \[\sqrt{n}(\hat{\theta}_{MLE}-\theta_0)\overset{D}{\rightarrow} N(0, I^{-1}_1(\theta_0)).\] When the dimension of the parameter is \(1\), a simple confidence interval is \[\hat{\theta}_{MLE} \pm z_{1-\alpha/2} \cdot \sigma_{\theta_0},\] where \(\sigma^2_{\theta_0} = I^{-1}_1(\theta_0)\) Such interval is not a confidence interval because \(\theta_0\) is unknown. We can modify it using a plug-in estimate of Fisher information: \[C_{n,\alpha} = [\hat{\theta}_{MLE} - z_{1-\alpha/2} \cdot \sigma_{\hat\theta_{MLE}}, \hat{\theta}_{MLE} + z_{1-\alpha/2} \cdot \sigma_{\hat\theta_{MLE}}],\] where \(z_\beta\) is the \(\beta\)-percentile of the standard normal distribution. Using Slutsky’s theorem, you can easily show that this confidence interval has asymptotic coverage.
When the dimension of the parameter is greater than \(1\), there are multiple ways we can construct a confidence interval. Note that in this case, the set \(C_{n,\alpha}\) is no longer an interval but a region/set so it is often called a confidence region/set. A simple approach of constructing a confidence set is via an ellipse. Note that asymptotic normality also implies (using the continuous mapping theorem) \[n (\hat{\theta}_{MLE}-\theta_0)^T I_1(\hat{\theta}_{MLE}) (\hat{\theta}_{MLE}-\theta_0)\overset{D}{\rightarrow} \chi^2_{p},\] where \(\chi^2_{p}\) denotes the \(\chi^2\) distribution with \(p\) degrees of freedom. So we construct the confidence set using \[C_{n,\alpha} = \left\{\theta: n (\hat{\theta}_{MLE}-\theta)^T I_1(\hat{\theta}_{MLE}) (\hat{\theta}_{MLE}-\theta)\leq \chi^2_{p,1-\alpha}\right\},\] where \(\chi^2_{p,\beta}\) is the \(\beta\)-percentile of the \(\chi^2\) distribution with \(p\) degrees of freedom.
Bootstrap confidence interval. The bootstrap approach is a Monte Carlo method for assessing the uncertainty of an estimator. It can be used to compute the variance of an estimator (not necessarily the MLE) and construct a confidence interval. In the case of likelihood inference, the bootstrap approach has the advantage that we do not need to know the closed form of \(I_1(\theta)\) to construct the confidence interval or to approximate the variance of the MLE.
While there are many variants of bootstrap methods, we introduce the simplest one – the empirical bootstrap. For simplicity, we assume that the dimension of \(\theta\) is \(1\) (the bootstrap works for higher dimensions as well). Let \(X_1,\cdots, X_n\) be the original sample. We then sample with replacement from the original sample to obtain a new sample of size \(n\), denoted \(X^*_1,\cdots, X^*_n\). This new sample is called a bootstrap sample. We find the MLE using the bootstrap sample and let \(\hat{\theta}_{MLE}^*\) denote the bootstrap MLE. Now we repeat the bootstrap process \(B\) times, leading to \(B\) bootstrap MLEs \[\hat{\theta}_{MLE}^{*(1)},\cdots, \hat{\theta}_{MLE}^{*(B)}.\] Let \(t_{\beta}\) denote the \(\beta\)-percentile of these \(B\) values, i.e., \[\hat {t}_\beta = \hat{G}^{-1}(\beta),\quad \hat{G}(t) = \frac{1}{B}\sum_{b=1}^B I(\hat{\theta}_{MLE}^{*(b)}\leq t).\] Then the bootstrap confidence interval of \(\theta_0\) is \[C_{n,\alpha} = [\hat{t}_{\alpha/2}, \hat{t}_{1-\alpha/2}].\] One can prove that under very mild conditions, the bootstrap confidence interval has asymptotic coverage.
The power of the bootstrap method is that we do not use anything about Fisher information! As long as we can compute the estimator, we can construct an asymptotically valid confidence interval. Note that if we do know the Fisher’s information, the bootstrap method can be modified using the bootstrap \(t\)-distribution method, which provides a better asymptotic coverage (namely, the \(o(1)\) decays faster to \(0\) than the above method and the normal confidence interval)2.
Statistical test is about how to design a procedure that allows us to make scientific discovery. Such a procedure has to be able to handle the uncertain nature of our data. In statistics, we model the data as random variables so the testing procedure needs to account for the randomness.
Let \(\mathcal{D}_n =\{X_1,\cdots, X_n\}\) denotes our data. The testing procedure involves two competing hypotheses:
Null hypothesis \(H_0\): the hypothesis that we want to challenge. It is often related to the current scientific knowledge.
Alternative hypothesis \(H_a\): the hypothesis that complements to the null hypothesis. It is the hypothesis we would like to prove to be plausible using our data.
The goal is to see if we have strong enough evidence (from \(\mathcal{D}_n\)) that we can argue the alternative hypothesis is more reasonable than the null hypothesis. If we do have enough evidence, then we will reject the null hypothesis. When the null hypothesis reflects the scenarios that can be explained by the current scientific knowledge, rejecting the null hypothesis means that we have discovered something new.
To design a testing procedure, we need to quantify the notion of evidence. The Frequentists and the Bayesian use different ways to measure the evidence. The Frequentist approach is the p-value whereas the Bayesian approach is the Bayes factor. We will talk about Bayes factor later so here we focus on the p-value.
Here is a summary on Frequentist approach of hypothesis test.
Based on the model and null hypothesis, design a test statistic.
Compute the distribution of the test statistics under the null hypothesis.
Plug-in the data into the test statistic, compute the probability of observing a more extreme data against the null hypothesis. Such a probability is the p-value.
Compare p-value to the significance level. If p-value is less than the significance level, we reject the null hypothesis.
The central idea of hypothesis test is to control the type-1 error, the probability of falsely rejecting \(H_0\) when \(H_0\) is correct. Essentially, the p-value can be interpreted as if we reject the null hypothesis (under this p-value), then our type-1 error is the same as the p-value. The significance level reflects the amount of type-1 error we can tolerate so when p-value is less than the significance level, we can reject \(H_0\). Due to the construction of p-value, a small p-value means that the null hypothesis does not fit to the data very well (so we are seeing an extreme event if \(H_0\) is true). Thus, small p-value or rejecting \(H_0\) under a small significance level means that we have more evidence against \(H_0\).
Note that there is another quantity called type-2 error, the probability of not rejecting \(H_0\) when \(H_0\) is false. Namely, type-2 error is concerned with the case that we fail to reject \(H_0\) when we should.
In statistics, we often control type-1 error first and the hope that the type-2 error is also small. When do we put more emphasis on type-1 error? This has something to do with the philosophy of scientific research. The scientific approach is a systematic way to acquire reliable knowledge. Thus, every discovery we made should be accompanied with sufficient evidences. In Frequentist approach, the measure of evidence against \(H_0\) is the p-value – the smaller p-value, the more evidence. Thus, controlling type-1 error means that we put requirements on the amount of evidence we need to claim a scientific discovery.
While there are many possible ways to construct a test statistic, here we consider two common approaches: Wald test and the likelihood ratio test.
Wald test. Assume in a simple case where we model the data as IID from a parametric model \(p(x;\theta)\) with a \(1D\) parameter. Suppose that we are comparing two hypotheses \[H_0: \theta = \theta_1,\qquad H_a: \theta\neq \theta_1.\] Since the MLE is a good estimate of \(\theta\), we can design our test statistic using the MLE. Because we know that the MLE has asymptotic normality, \[\sqrt{n}I_1(\theta)^{-1/2} (\hat{\theta}_n-\theta)\overset{D}{\rightarrow} N(0,1),\] we can then use \[T_n(\theta) = \sqrt{n}I_1(\theta)^{-1/2} (\hat{\theta}_n-\theta)\] as our test statistic. If \(H_0\) is true, \(T_n(\theta_1)\) should behaves like a standard normal distribution. When we observe the actual data, we can then compare the observed value of \(T_n(\theta_1)\) against \(0\) (since \(T_n(\theta_1) = 0\) means a perfect match with \(H_0\)). We reject \(H_0\) if \(T_n(\theta_1)\) is either too large or too small. The p-value is then \[p(\mathcal{D}_n) = 1-2\cdot \Phi^{-1}(|T_n(\theta_1)|),\] where \(\Phi(t) = P(Z\leq t)\) is the CDF of a standard normal distribution. Note that the above method is also known as the Wald test.
Likelihood ratio test (LRT). LRT is another popular way to conduct hypothesis tests under a parametric model. It can be easily applied to multivariate parameters so now we assume that there are \(p\) parameters, i.e., \(\theta\in \Theta\subset \mathbb{R}^p\). The null hypothesis and the alternative hypothesis are \[H_0: \theta\in\Theta_0,\qquad H_a: \theta\in \Theta\backslash \Theta_0.\] Let \(\ell_n(\theta)\) be the log-likelihood function and let \[\hat\theta_{MLE} = {\sf argmax}_{\theta\in\Theta} \ell_n(\theta),\quad \hat{\theta}_{MLE,0} = {\sf argmax}_{\theta\in\Theta_0}\ell_n(\theta)\] be the global MLE and the constrained MLE (in the null hypothesis space). The LRT starts with a test statistic \[{\sf LRT} = \frac{\sup_{\theta\in\Theta_0} L_n(\theta)}{\sup_{\theta\in\Theta}L_n(\theta)} = \frac{L_n(\hat{\theta}_{MLE,0})}{L_n(\hat{\theta}_{MLE})}\] and rejects \(H_0\) if \({\sf LRT}\) is too small. To get the p-value of the LRT, we often use the following fact: \[-2\log {\sf LRT}= 2\left(\ell_n(\hat\theta_{MLE} ) - \ell_n(\hat\theta_{MLE,0} )\right)\overset{D}{\rightarrow} \chi^2_{r},\] where \(r = {\sf dim}(\Theta)-{\sf dim}(\Theta_0)\) is the difference in dimensions (degrees of freedom) between the full parameter space and constrained parameter space. The above result is also known as the likelihood ratio approximation or Wilks’ theorem. The LRT proposes to use \(T_n = 2\left(\ell_n(\hat\theta_{MLE} ) - \ell_n(\hat\theta_{MLE,0} )\right)\) as our test statistic and compare it with the CDF of \(\chi^2_{r}\) to obtain the p-value.
Remarks.
Why \(r\) degrees of freedom? Some of you may be wondering why we are getting a \(\chi^2\) distribution with \(r\) degrees of freedom. Here is a simple explanation using geometry. Recall that the MLE behaves like a normal distribution around the true parameter \(\theta_0\). If \(H_0\) is true, \(\theta_0\in\Theta_0\) so \(\hat\theta_{MLE}\) will be asymptotically normally distributed around \(\theta_0\). Since \(\hat\theta_{MLE}\) is the maximizer over \(\Theta\), it has \(p\) degrees of freedom (it can move in each of the \(p\) dimensions). The constrained MLE \(\hat\theta_{MLE,0}\) is the maximizer under \(\Theta_0\), which has \(r\) constraints. Thus, \(\Theta_0\) has \(p-r\) degrees of freedom, which implies that its maximizer \(\hat\theta_{MLE,0}\) uses \(p-r\) degrees of freedom. The remaining number of degrees of freedom is \(p-(p-r) = r\) so this is why we obtain a \(\chi^2\) distribution with \(r\) degrees of freedom. Note that one can replace the \(r\) constraints by saying that \(\Theta_0\) is a \((p-r)\)-dimensional manifold.
Equivalence between likelihood ratio test and Wald test. Asymptotically, one can show that the likelihood ratio test and the Wald test are equivalent under classical assumptions on the MLE. Also, there is another test that is closely related to them called the score test, which is based on the value of the score function as a test statistic. Again, asymptotically the score test and the other two tests are equivalent. However, when the likelihood function is more complex, such as having multiple local maxima, these three tests may not be the same (they can be quite different)3.
Relation to confidence interval. Hypothesis test can be used to construct a confidence interval. Consider testing the null hypothesis: \(H_0:\theta = \theta_1\) for a specific value \(\theta_1\) under a significance level \(\alpha\). For each \(\theta_1\in\Theta\), we can do a hypothesis test. Some parameters will be rejected whereas the others will not. Let \(\hat{A}_{n,\alpha}\) be the collection of parameters that the null hypothesis will not be rejected. Then you can show that \(\hat{A}_{n,\alpha}\) is a confidence interval of \(\theta\). This approach is also known as the confidence interval from inverting a hypothesis test. Note that if we are given a \(1-\alpha\) confidence interval of the parameter of interest \(\theta\), we can use it to test the null hypothesis \(H_0: \theta\in \Theta_0\). If the confidence interval intersects with \(\Theta_0\), we cannot reject \(H_0\). If they are disjoint, then we can reject \(H_0\) under a significance level \(\alpha\). Thus, hypothesis test problem and confidence intervals are highly related to each other. However, testing a given hypothesis is often an easier problem than constructing a confidence interval because confidence interval requires a procedure that is valid regardless of the null hypothesis being correct or not. On the other hand, hypothesis test problem only requires a procedure that works when \(H_0\) is true.
There are many misconceptions about p-values. To clarify what p-value stands for and what it does NOT stand for, I obtain the following 6 principles from American Statistical Association’s website4:
P-values can indicate how incompatible the data are with a specified statistical model.
P-values do not measure the probability that the studied hypothesis is true, or the probability that the data were produced by random chance alone.
Scientific conclusions and business or policy decisions should not be based only on whether a p-value passes a specific threshold.
Proper inference requires full reporting and transparency.
A p-value, or statistical significance, does not measure the size of an effect or the importance of a result.
By itself, a p-value does not provide a good measure of evidence regarding a model or hypothesis.
Many theory about the MLE assumes that the population distribution function belongs to our parametric family. However, this is a very strong assumption in reality. It is very likely that the population distribution function does not belong to our parametric family (e.g., the population PDF is not Gaussian but we fit a Gaussian to it). What will happen in this case for our MLE? Will it still converge to something? If so, what will be the quantity that it is converging?
Model mis-specification studies the situation like this – we assume a wrong model for the population distribution function. Let \(p_0(x)\) be the population PDF and we assume that the population PDF can be written as \(p(x;\theta)\). However, \(p_0\neq p(x;\theta)\) for every \(\theta\in \Theta\). It turns out that the MLE \(\hat{\theta}_{MLE}\) still converges under mild assumptions to a quantity \(\theta^*\) in probability. Moreover, the corresponding PDF/PMF \(p(x;\theta^*)\) has an interesting relation with \(p_0(x)\). Assume that the RV \(X\) has a PDF/PMF \(p_0\). Then \[\mathbb{E}\left\{\log \left(\frac{p_0(X)}{p(X;\theta^*)}\right)\right\} = \inf_{\theta\in\Theta} \mathbb{E}\left\{\log \left(\frac{p_0(X)}{p(X;\theta)}\right)\right\}= \inf_{\theta\in\Theta} {\sf KL}(p_0, p_\theta),\] where \({\sf KL}\) is also known as the Kullback-Liebler (KL) divergence and \(p_\theta(x) = p(x;\theta)\). Namely, the MLE corresponds to the parametric distribution in the specified family that minimizes the KL divergence to the population distribution.
In the model mis-specification case, the MLE still satisfies the score equation (under appropriate assumptions) but the Fisher information may not reflect the actual curvature of the likelihood function around \(\theta^*\). The asymptotic covariance (related to the curvature) of \(\hat{\theta}_{MLE}\) will be \(\Sigma = I^{-1}_1(\theta^*) \mathbb{E}(S_1(\theta^*)S^T_1(\theta^*))I^{-1}_1(\theta^*)\) and we still have \[\sqrt{n}(\hat{\theta}_{MLE} - \theta^*) \overset{D}{\rightarrow}N(0,\Sigma).\]
The Bayesian inference is an alternative statistical paradigm to the Frequentist approach. The Bayesian approach interprets the probability in a broader sense that include subjective probability, which allows us to assign probability to almost every quantity in our model (including the parameter of interest and even a statistical model). The Bayesian inference relies on a simple decision theoretic rule – if we are competing two or more choices, we always choose the one with higher probability. This simple rule allows us to design an estimator, construct an interval, and perform hypothesis test.
In Bayesian analysis, we assign a probability to every parameter in our model. For a parametric model \(p(x;\theta)\), the parameter of interest \(\theta\) is given a prior distribution \(\pi(\theta)\) that reflects our belief about the value of \(\theta\). In a sense, the prior distribution quantifies our subjective belief about the parameter \(\theta\). A higher value of \(\pi(\theta)\) indicates that we believe \(\theta\) is a more likely value.
How do we interpret this prior distribution? Here is a decision-theoretic way of viewing it. To simplify the problem we assume that \(\Theta = \{0,1,2\}\). Even without any data at hand, we can ask ourselves about our belief about each parameter value. Some people may think that \(1\) is the most likely one; some may think that \(2\) is the most likely one. To make our belief more precise, we use probability to quantify it. Let \(\pi(j)\) be the number that reflects our belief about \(\theta = j\). We interpret the numerical value of \(\pi(j)\) as follows. We are forced to guess between \(\theta=j\) versus \(\theta\neq j\). If the answer is \(\theta=j\) and we indeed guess correctly, we will be rewarded \(\delta\) dollars. If the answer is \(\theta\neq j\) and we get it correct, we will be rewarded \(1\) dollar. If we get it wrong, we do not lose anything. Our principle is to maximize our expected reward. Now assume that the true value of \(\theta\) has equal probability of being \(j\) or not \(j\). Then what should we choose? \(\theta=j\) or \(\theta\neq j\)? Now we think about this problem by varying \(\delta\) from \(0\) to infinity. When \(\delta\) is small, unless we have a very strong belief in \(\theta=j\), we will not bid on it. When increasing \(\delta\), at a certain threshold we will switch our decision from bidding on \(\theta\neq j\) to \(\theta=j\). Let this threshold be \(\eta_j\). \(\eta_j\) is a number that reflects our belief about \(\theta = j\) and we associate it with our prior \[\pi(j) = \frac{1}{1+\eta_j}\Leftrightarrow \eta_j = \frac{1-\pi(j)}{\pi(j)}\quad\mbox{(odds of $\theta=j$)}.\] Here, you see that we only use one simple decision rule – bidding on the choice with a higher expected outcome. This allows us to quantify our belief.
Using the prior distribution, the Bayesian probability model can be written as follows: \[\begin{align*} X_1,\cdots, X_n|\theta&\overset{IID}{\sim} p(x|\theta)\\ \theta&\sim \pi. \end{align*}\] The Bayesian inference focuses on the distribution of \(\theta\) after observing \(X_1,\cdots,X_n\): \[\pi(\theta|X_1,\cdots,X_n) = \frac{p(X_1,\cdots, X_n,\theta)}{p(X_1,\cdots, X_n)} \propto \underbrace{p(X_1,\cdots, X_n|\theta)}_\text{likelihood} \times \underbrace{\pi(\theta)}_\text{prior}.\] This distribution is also known as the posterior distribution.
The posterior distribution informs us about how our prior belief is updated after seeing the data. It is the central quantity in Bayesian inference – all our decisions will be related to it. From a Bayesian’s point of view, probability models are just mathematical tools for analyzing data. We do not assume that data are generated from a true probability distribution. We just view the data as generated from \(p(x;\theta)\). Given that we do not assume the probability model to be the true model, there is NO true parameter so we cannot talk about conventional statistical errors. However, Bayesians do have another way of expressing uncertainty in our inference – the posterior distribution. The posterior distribution reflects our belief about the parameter after seeing the data, we can use it as a measure of uncertainty about \(\theta\). If the posterior distribution is more spread out, then the uncertainty in our inference is larger. On the other hand, if the posterior distribution is very concentrated, then there is very little (Bayesian) uncertainty.
There are two common estimators in Bayesian inference: the posterior mean and maximum a posteriori estimation (MAP).
Posterior mean. Just like we often use the sample mean as an estimator of the population mean, the mean of the posterior distribution is a common quantity that is used as an estimator of \(\theta\): \[\hat{\theta}_\pi = \mathbb{E}(\theta|X_1,\cdots, X_n)= \int \theta \cdot\pi(\theta|X_1,\cdots,X_n) d\theta.\] It represents the average location of our belief about the parameter after seeing the data.
Maximum a posteriori estimation (MAP). Another common estimator of \(\theta\) is the MAP; it relies on a similar principle as the MLE – we choose the one that is the most likely. Here ’likely’ is interpreted as our posterior belief about the parameter of interest \(\theta\). Formally, MAP is defined as \[\hat{\theta}_{MAP} = {\sf argmax}_\theta \pi(\theta|X_1,\cdots, X_n).\]
Example: Binomial Sampling. Assume that we have an observation \(Y\sim {\sf Bin}(N,\theta)\) where \(N\) is known and the parameter of interest is \(\theta\): \[P(Y=y|\theta) = {N\choose y} \theta^y(1-\theta)^{N-y}.\] We use a Beta distribution with parameters \((\alpha,\beta)\) as our prior distribution for \(\theta\). Namely, \[\pi(\theta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^{\alpha-1}(1-\theta)^{\beta-1},\] where \(\Gamma(z) = \int_0^\infty t^{z-1}e^{-t}dt\) is the Gamma function and \(\alpha,\beta>0\). Note that \((\alpha,\beta)\) are called the hyperparameters and are known quantities (because we know our belief about the data). For a Beta distribution with parameters \(\alpha,\beta\), the mean is \(\frac{\alpha}{\alpha+\beta}\).
The posterior distribution is \[\begin{align*} \pi(\theta|Y) &= \frac{{N\choose Y} \theta^Y(1-\theta)^{N-Y} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^{\alpha-1}(1-\theta)^{\beta-1}}{ \int {N\choose Y} \theta^Y(1-\theta)^{N-Y} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^{\alpha-1}(1-\theta)^{\beta-1} d\theta}\\ & \propto\theta^{Y+\alpha-1}(1-\theta)^{N-Y+\beta-1} \end{align*}\] so it is a Beta distribution with parameters \((Y+\alpha, N-Y+\beta)\). Then the posterior mean and MAP are \[\hat{\theta}_\pi = \frac{Y+\alpha}{N+\alpha+\beta},\qquad \hat{\theta}_{MAP} = \frac{Y+\alpha-1}{N+\alpha+\beta-2}\] (these are the mean and the mode of a Beta distribution).
Note that in this problem, the MLE is \(\hat{\theta}_{MLE} = \frac{Y}{N}\). Thus, the posterior mean has an interesting decomposition: \[\begin{align*} \hat{\theta}_\pi &= \frac{Y+\alpha}{N+\alpha+\beta} \\ &= \hat{\theta}_\pi = \frac{Y}{N+\alpha+\beta} +\frac{\alpha}{N+\alpha+\beta} \\ & = \frac{Y}{N} \times \frac{N}{N+\alpha+\beta}+ \frac{\alpha}{\alpha+\beta} \times \frac{\alpha+\beta}{N+\alpha+\beta}\\ & = \hat{\theta}_{MLE} \times W + \mbox{Prior mean} \times (1-W), \end{align*}\] where \(W = \frac{N}{N+\alpha+\beta}\) is a weight that is tending to \(1\) when \(N\rightarrow \infty\). This phenomenon – the posterior mean can be written as the weighted average of the MLE and the prior mean – occurs in several scenarios. Moreover, the fact that the weights \(W\rightarrow 1\) as the sample size \(N\rightarrow\infty\) means that when we have more and more data, the prior distribution seems to be irrelevant. Thus, the posterior mean would have a similar asymptotic property as the sample mean. However, this is not a general phenomenon; often only certain combination of prior and likelihood models will have this feature.
Remarks.
Choice of prior and conjugate prior. The choice of prior reflects our belief about the parameter before seeing any data. Sometimes people want to choose a prior distribution such that the posterior distribution is in the same family as the prior distribution, just like what we have observed in the above example. If a prior distribution and a likelihood function lead to a posterior that belongs to the same family as the prior, we call this prior a conjugate prior. There are several conjugate priors known to date, see https://en.wikipedia.org/wiki/Conjugate_prior for an incomplete list.
Another common choice of prior is called the Jeffreys prior5, which chooses a prior \(\pi(\theta)\propto \sqrt{{\sf det}(I_1(\theta))}\), where \(I_1(\theta)\) is the Fisher information matrix. One can view the Jeffreys prior as the prior that we do not have any prior belief about \(\theta\); or more formally, an uninformative prior.
Challenge of computing the posterior. In general, if we do not choose a conjugate prior, the posterior distribution could be difficult to compute. The challenge often comes from the normalization quantity \(p(X_1,\cdots, X_n)\) in the denominator of the posterior \(\pi(\theta|X_1,\cdots,X_n)\) (the numerator is just the prior times the likelihood). In practice we will use Monte Carlo method to compute the posterior – we generate points from \(\pi(\theta|X_1,\cdots,X_n)\) and as we generate enough points, these points should approximate the true posterior distribution well. We will talk more about this later in the lecture of MCMC (Monte Carlo Markov Chain).
Consistency. In pure Bayesian’s point of view, statistical consistency is not an important property because probability model is a working model to describe the data and we do not need to assume that there exists an actual parameter that generates the data. Thus, the posterior distribution is the quantity that we really need to make our inference. However, sometimes Bayesian estimators, such as the posterior mean or MAP, does have statistical consistency. Namely, \(\hat{\theta}_\pi \overset{P}{\rightarrow}\theta_0\) and \(\hat{\theta}_{MAP}\overset{P}{\rightarrow}\theta_0\), where the data \(X_1,\cdots,X_n\overset{IID}{\sim} p(x;\theta_0).\) This is often related to the Bernstein-von Mises theorem6. Although statistical consistency was not an important property in Bayesian paradigm (because Bayesian does not assume the data is indeed from a probability model; probability models are just a mathematical model to help us analyze the data), still many researchers would prove consistency when proposing a Bayesian approach.
The Bayesian interval is very straightforward – since the posterior contains our belief about the parameter after seeing the data, we just construct the interval using the posterior distribution. Given a credible level \(1-\alpha\), a credible interval \(D_{n,\alpha}\) satisfies \[1-\alpha = \int_{D_{n,\alpha}}\pi(\theta|X_1,\cdots, X_n)d\theta.\] Often we will choose the credible interval such that it is the shortest. It turns out that (actually it is not hard to prove) such a credible interval is related to the upper level set of the posterior. Given a level \(\lambda\), we define \[L_\lambda = \{\theta\in \Theta: \pi(\theta|X_1,\cdots,X_n)\geq \lambda\}.\] Define \(V(\lambda) = \int_{L_\lambda}\pi(\theta|X_1,\cdots, X_n)d\theta\). When \(\lambda\) is very large, \(V(\lambda)\) is very small. When we decrease \(\lambda\), \(V(\lambda)\) will increase since we are including more regions. At a critical level \(\lambda_\alpha\), we will have exactly \[V(\lambda_\alpha) = 1-\alpha.\] We will then use the level set \(L_{\lambda_\alpha}\) as the credible interval. You can show that under appropriate conditions (\(\pi(\theta|X_1,\cdots,X_n)\) has no flat regions and is Lipschitz), \(L_{\lambda_\alpha}\) is the shortest credible interval with a credible level \(1-\alpha\).
Because the above (shortest) credible interval is very straightforward, when people talk about credible intervals, they are often referring to this interval.
Bayesian hypothesis testing is also very straightforward – as you can guess, it is based on the posterior distribution. Instead of just putting priors on parameters, using our decision-theoretic way, we can also put priors on the hypotheses.
Recall that \(H_0\) and \(H_a\) are the null and alternative hypotheses. Let \(\pi(H_0)\) and \(\pi(H_a) = 1-\pi(H_0)\) denote the prior distribution on the two hypotheses and let \(\pi(H_0|X_1,\cdots,X_n)\) and \(\pi(H_a|X_1,\cdots, X_n)\) be the posterior distribution given the data.
The testing procedure is very simple – we reject \(H_0\) if \[\begin{equation} \frac{\pi(H_0|X_1,\cdots,X_n)}{\pi(H_a|X_1,\cdots,X_n)} <1\Leftrightarrow \frac{p(X_1,\cdots,X_n|H_0)}{p(X_1,\cdots,X_n|H_a)} \cdot \frac{\pi(H_0)}{\pi(H_a)} <1. \end{equation}\] Namely, we reject \(H_0\) if our posterior belief about \(H_0\) is less than that of \(H_a\).
In many scenarios, the hypotheses will not directly give us a probability model related to the data. In a parametric model, they often put some constraints on the parameter. Thus, we can then rewrite the posterior \(p(X_1,\cdots,X_n|H_0)\) as \[p(X_1,\cdots,X_n|H_0) = \int p(X_1,\cdots,X_n, \theta|H_0)d\theta = \int p(X_1,\cdots,X_n|\theta,H_0) \pi(\theta|H_0) d\theta.\] Using the above equality, we define the Bayes factor as \[\begin{equation} {\sf BF} (X_1,\cdots, X_n) = \frac{p(X_1,\cdots,X_n|H_0)}{p(X_1,\cdots,X_n|H_a)} = \frac{\int p(X_1,\cdots,X_n|\theta,H_0) \pi(\theta|H_0) d\theta}{\int p(X_1,\cdots,X_n|\theta,H_a) \pi(\theta|H_a) d\theta}. \end{equation}\] Using the Bayes factor, we can rewrite the decision rule as: we reject \(H_0\) if \[{\sf BF} (X_1,\cdots, X_n) \cdot \underbrace{\frac{\pi(H_0)}{\pi(H_a)}}_\text{odds} <1.\] Namely, if the Bayes factor is less than the inverse of the odds of \(H_0\), we reject \(H_0\).
The Bayes factor can be viewed as a Bayesian version of p-value – the smaller, the less favor for \(H_0\). To see how the Bayes factor is like the p-value, note that the Frequentist way of rejecting \(H_0\) is if the p-value is less than a pre-specified significance level \(\alpha\). The Bayesian’s threshold is given by the odds of \(H_0\) from the prior distribution. If the inverse of the Bayes factor is greater than the odds of \(H_0\), we reject \(H_0\). In Bayesian, the threshold in testing (significance level) has a simple interpretation – the odds of our prior belief about the null hypothesis.
Example: Binomial Sampling. Now we come back to the Binomial sampling example where we have an observation \(Y\sim {\sf Bin}(N,\theta)\) where \(N\) is known and the parameter of interest is \(\theta\): \[P(Y=y|\theta) = {N\choose y} \theta^y(1-\theta)^{N-y}.\] We use a Beta distribution with parameters \((\alpha,\beta)\) as our prior distribution for \(\theta\). Now we consider testing \[H_0: \theta = 0.5,\qquad H_a:\theta\neq 0.5.\] Our goal is to compute the Bayes factor. First, we compute the numerator: \[p(Y|H_0) = \int_0^1 p(Y|\theta,H_0) \pi(\theta|H_0) d\theta = \int_0^1p(Y|\theta) \delta(\theta= 0.5) d\theta = p(Y|\theta=0.5) = {N\choose Y} 0.5^Y0.5^{N-Y} = {N\choose Y} 0.5^N.\] The denominator of the Bayes factor will be \[\begin{align*} p(Y|H_a) &=\int_0^1 p(Y|\theta,H_a) \pi(\theta|H_a) d\theta\\ & = \int_0^1 p(Y|\theta,H_a) \pi(\theta) d\theta \qquad \mbox{(Note that $\theta=0.5$ is just a single point so it does not affect the integral)}\\ & = \int_0^1{N\choose Y} \theta^Y\theta^{N-Y} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^{\alpha-1}(1-\theta)^{\beta-1}d\theta \\ & = {N\choose Y} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma(Y+a)\Gamma(N-Y+b)}{\Gamma(N+a+b)}. \end{align*}\] Therefore, the Bayes factor is \[{\sf BF}(Y) = \frac{0.5^N}{\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma(Y+a)\Gamma(N-Y+b)}{\Gamma(N+a+b)}}.\]
As you can see, the computation of Bayes factor does not involve the prior on \(H_0\). However, it does involve the prior of the parameter \(\theta\) (and this is why the value of it depends on the hyperparameters \(\alpha,\beta\)).
There is an interesting relation between the Bayes factor and the likelihood ratio test. In the above Binomial example, if we use the likelihood ratio test, we obtain a test statistic: \[{\sf LRT} = \frac{p(Y|\theta=0.5)}{p(Y|\theta=\hat{\theta}_{MLE})} = \frac{{N\choose Y} 0.5^N}{\sup_{\theta\in[0,1]} {N\choose Y} \theta^Y(1-\theta)^{N-Y}} = \frac{0.5^N}{ (Y/N)^Y(1-Y/N)^{N-Y}}.\] The numerator is the same while the denominator is slightly different – the LRT only uses the maximum likelihood value whereas the Bayes factor uses the average value within \(H_a\)! An interesting fact: if the alternative hypothesis is a simple hypothesis \(H_a: \theta=\theta_a\) for a fixed quantity \(\theta_a\), then the Bayes factor and the \({\sf LRT}\) statistic coincide.
The Bayesian approach can also be applied to nonparametric estimation. Nonparametric estimation refers to the case where we do not assume the parameter of interest is a vector or number. One common problem is density estimation. If our goal is to estimate the underlying PDF or CDF, we can do it without assuming the data are from a parametric family such as a Gaussian. For instance, a histogram can be used as a density estimator without assuming the data is from a Gaussian or any other parametric family (there are many more advanced techniques such as the kernel density estimator, orthonormal basis approach, wavelet approach, ...etc).
The challenge of Bayesian inference in nonparametric problems is that we need to put a prior on a ‘function space’ (a collection of functions). Often a function space does not admit a density function. However, there are some tricks so that we can put a prior on a function space.
Dirichlet process. The Dirichlet process7 is a stochastic process that generates random probability distributions. It also has several other interesting names such as the Chinese buffet process and stick-breaking process. The appealing feature of the Dirichlet process is that it is constructed from an algorithmic scheme so there is a simple way to sample points from a random distribution generated from this process.
Sequential mixture model. Another way to assign priors on distributions is to use a sequential mixture model, for instance, a Gaussian mixture model with parameters drawn from some distribution. When we allow the number of mixture components to increase with sample size, we can approximate many distributions in the function space. Note that when we allow the number of mixture to increase, this implies that our prior distribution is changing with respect to the sample size.
Basis approach. A smooth density function with a compact support may be written as \(p_0(x) = \sum_{k=1}^\infty \theta_k \phi_k(x)\), where \(\{\phi_k(x): k=1,2,\cdots\}\) is a basis and \(\{\theta_k: k=1,2,\cdots\}\) are the coefficients. For instance, one may choose the cosine basis or spline basis. Putting priors on \(\theta_k\) leads to a prior distribution over functions. Note that in practice we often had to truncate the basis at certain level , i.e., we use only \(\{\phi_k(x): k=1,2,\cdots, N\}\) for some \(N=N\). Often we allow \(N= N_n\) to increase with respect to the sample size so that our prior covers most part of the function space.
\(\epsilon\)-cover/bracketing. There is another approach to assign priors over function space using the covering/bracketing of function space. An \(\epsilon\)-cover is a collection of functions in a function space such that every function inside the space has a distance at most \(\epsilon\) to the nearest element inside the cover (bracketing is a generalization of this concept). There are many \(\epsilon\)-covers but we often use those with the minimal number of elements. Now consider a sequence of \(\epsilon\)-cover with \(\epsilon=\epsilon_n\rightarrow 0\) as \(n\rightarrow\infty\). Let \(\Pi_n\) be a uniform distribution over each element of \(\epsilon_n\)-cover. Then the mixture distribution \(\Pi = \sum_{n=1}^\infty \beta_n\Pi_n\), with an appropriate choice of \(\{\beta_n\}\) forms a prior distribution that can well-approximate almost every function in the function space (when the function space is sufficiently smooth). This idea is a powerful tool in constructing a prior with amazing theoretical properties8 although it is often hard to numerically compute this prior.
Prediction is about the question that given observations \(X_1,\cdots,X_n,\) what will be the possible range of the next observation \(X_{n+1}\)? Using probability model, we know that the next observation \(X_{n+1}\) will have a distribution. So the prediction problem is about finding the distribution \(p(x_{n+1}|X_1,\cdots,X_n)\). Here we will show how the Frequentist prediction method and the Bayesian prediction method are different in a parametric model.
The Frequentist approach is very simple – given that we have already assigned a parametric model \(p(x;\theta)\) for the observations, we should just use it as our predictive distribution. To use this model, we need to choose the parameter \(\theta\). A simple choice is the MLE \(\hat{\theta}_{MLE}\). Thus, we will predict the distribution of \(X_{n+1}\) as \[p(x_{n+1}|X_1,\cdots,X_n) = p(x_{n+1};\theta= \hat{\theta}_{MLE}).\]
The Bayesian approach again relies on the posterior. The predictive distribution \(p(x_{n+1}|X_1,\cdots,X_n)\) can be written as \[\begin{align*} p(x_{n+1}|X_1,\cdots,X_n) & =\int p(x_{n+1},\theta|X_1,\cdots,X_n)d\theta\\ & = \int p(x_{n+1}|\theta,X_1,\cdots,X_n) \pi(\theta|X_1,\cdots,X_n)d\theta\\ & = \int p(x_{n+1}|\theta) \underbrace{\pi(\theta|X_1,\cdots,X_n)}_\text{posterior distribution}d\theta. \end{align*}\] Thus, the predictive distribution is the averaged distribution of \(p(x|\theta)\) where we average \(\theta\) over the posterior distribution.
Here, as you can see, the two paradigms make prediction using different principles – the Frequentists use only the most likely model to make predictions whereas the Bayesians use the averaged model over the posterior distribution to make predictions.
Both Frequentist and Bayesian approaches are self-consistent. They start with probability models and design their own procedure for estimation, interval inference, and hypothesis test. In practice, it is hard to really say if any method truly describes the reality because we do not even know if our data are indeed generated from a random process. Some people believe that our data should be viewed as realizations from a complex dynamic system and they can still construct estimators and derive consistency without introducing any probability (often they would use the Ergodic theory9). Thus, we cannot say if Frequentist or Bayesian approach is the right approach to analyze data – they are just principles that allows us to analyze data in a well-established way. In what follows, I briefly comments on the criticisms and defends of the two paradigms.
Many people support Frequentist paradigm because it is more objective – we do not introduce any subjective belief (prior) on the parameter of interest. Moreover, the way Frequentist views the parameter of interest, an unknown but fixed quantity, fits into the what most scientists think about parameters in model – these parameters are some fixed numbers but we just do not know it. Moreover, the Frequentist’s view of the probability – long term frequency – is very intuitive for most people.
The Bayesian paradigm is very clean. From estimation, interval inference, to hypothesis test, every inference just depends on one single quantity – the posterior distribution. And we only have one single guiding principle – decision theoretic rule – we choose the one that has a higher posterior distribution. On the other hand, in the case of Frequentist approach, estimation requires a principle (such as the MLE principle), confidence interval relies on another principle (coverage), and the hypothesis test introduces another principle (p-value and significance level). So Bayesian is a very elegant way to analyzing the data.
A major criticism of Bayesian is the concept of subjectivity – the prior distribution. Many people prefer Frequentist approaches over Bayesian approaches because they think that scientific studies should be objective. However, this argument is actually not valid if we really think deep about Frequentist approach – when we use a probability model to describe the data, we are already making a subjective choice of how we model the data! Why not use a dynamic system approach? why not use a Bayesian? The choice of Frequentist approach itself is a subjective decision made by scientists. Moreover, the choice of estimator, the choice of confidence level, and the choice of significance level along with the testing procedures, are all subjective decisions. None of them are truly objective. One attractive feature of Bayesian paradigm is that Bayesians not only accept the fact that we are making many subjective choices in analyzing data but also they have a well-defined mathematical framework – the probability model – to describe how these subjectivity are propagating throughout the analysis.
In both the likelihood model and the Bayesian approach, our estimator is often obtained by a maximization/minimization procedure such as the MLE. In simple problems, the MLE has a closed-form solution. But in many practical scenarios, the MLE does not have a closed form and has to be computed via a numerical procedure. In this section, we discuss a generic but popular approach for numerically computing the MLE: the gradient descent method.
Gradient descent is a generic approach for finding a minimum of a smooth function. Note that finding the minimum and finding the maximum are almost equivalent questions (just flip \(f\) to \(-f\) or \(1/f\)) so here we will focus on finding the minimum. Let \(f\) be a smooth function. Our goal is to find \(x^* = {\sf argmin}_x f(x)\), the location where the minimum of \(f\) occurs, and the minimal value \(f^* = \min_x f(x)\).
In gradient descent (and almost all other minimization methods), the shape of a function is very important. When the function has a nice shape (like a Gaussian), there is a unique and well-defined minimum and gradient descent is able to find it in a reasonable time.
We will start with the simplest case by assuming that \(f\) is a smooth and convex function. Recall that convexity implies that for any \(x,y\) and \(\alpha\in [0,1]\), \[f(\alpha x+(1-\alpha) y) \leq \alpha f(x ) + (1-\alpha) f(y).\] The gradient descent algorithm works as follows.
Choose an initial point \(x^{(0)}\) and a step size (learning rate) \(\gamma>0\).
Iterates the following until some convergence criterion is met: \[x^{(t+1)} \leftarrow x^{(t)} - \gamma \nabla f(x^{(t)}).\]
Namely, each time we move from the current point \(x^{(t)}\) by \(\gamma\) times the value of gradient. The following picture shows a gradient descent path starting from 4 different initial points
The step size \(\gamma\) is very important in gradient descent. It affects the convergence in two senses. If \(\gamma\) is too large, the gradient descent may diverge. If \(\gamma\) is too small, the gradient descent may converge very slowly. Note that one can also choose \(\gamma = \gamma_t\) that depends on the number of steps.
In the top left panel, we choose a very large step size and the gradient descent diverges at the end. In the top right panel, we choose a very small step size and the algorithm converges very slowly. Therefore, how to choose a good step size is very important. There are two common strategies for this.
Backtracking line search.
Choose \(0<\beta<1\) and \(0<\alpha\leq 1/2\).
At each iteration, start with \(\gamma = \gamma_0\), and while \[f(x^{(t)}-\gamma\nabla f(x^{(t)})) > f(x^{(t)}) - \alpha \gamma\|\nabla f(x^{(t)})\|^2_2,\] shrink \(\gamma = \beta \gamma\). Else, perform gradient descent update \(x^{(t+1)} \leftarrow x^{(t)} - \gamma\nabla f(x^{(t)})\).
Note that often we will choose \(\alpha=1/2\).
Exact line search. The exact line search attempts to choose \(\gamma\) based on \[\gamma = {\sf argmin}_{s\geq 0} f(x^{(t)}-s\nabla f(x^{(t)})).\] It is another optimization problem so it is hard to obtain the minimizer exactly (but we can do approximation to it).
Note that in both backtracking and exact line search, one has to choose the step size at every time step \(t\).
Here our convergence analysis is not about the statistical convergence but the algorithmic convergence. We want to analyze how close is our minimizer \(x^{(t)}\) to the true minimum \(x^*\). We will analyze the convergence rate in terms of \(t\), the number of steps. Let \(f^* = f(x^*) = \min_x f(x)\).
Before stating our main result, we first note the following two important properties of \(f\) being \(L\)-smooth and convex. \[\begin{align} \mbox{($L$-smooth) }\|\nabla f(x) - \nabla f(y)\|_2 \leq L\|x-y\|_2&\Rightarrow f(x) \leq f(y) + (x-y)^T \nabla f(y)+ \frac{L}{2}\|y-x\|^2 \label{eq::Lsmooth}\\ \mbox{(convexity) }f\mbox{ is convex }&\Rightarrow f(x) \leq f(y) + (x-y)^T \nabla f(x) \label{eq::convexity} \end{align}\] Note that the property \[\|\nabla f(x) - \nabla f(y)\|_2 \leq L\|x-y\|_2\] is also known as \(L\)-smooth.
For a function that is \(L\)-smooth and convex, we have the following useful lemma.
Lemma 1. Suppose that \(f\) is \(L\)-smooth and convex and the step size \(\gamma\leq 2/L\), then \[\begin{align} f(y)-f(x)&\leq (y-x)^T\nabla f(y) - \frac{1}{2L}\|\nabla f(y)-\nabla f(x)\|^2 \label{eq::lem1}\\ (y-x)^T (\nabla f(y)-\nabla f(x)) &\geq \frac{1}{L}\|\nabla f(x)-\nabla f(y)\|^2 \label{eq::lem2}\\ f(x^{(t+1)})&\leq f(x^{(t)})- \frac{\gamma}{2}(2-\gamma L)\|\nabla f(x^{(t)})\|^2 \label{eq::lem3}\\ \|x^{(t+1)}-x^*\|^2&\leq \|x^{(t)}-x^*\|^2- \frac{\gamma}{L}(2-\gamma L)\|\nabla f(x^{(t)})\|^2 \label{eq::lem4}\\ f(x^{(t)})-f^* &\leq \|x^{(0)}-x^*\| \|\nabla f(x^{(t)})\| \label{eq::lem5} \end{align}\]
With the above lemma, we can prove the following convergence rate of a convex and \(L\)-smooth function.
Theorem 2. Assume that \(f\) is convex and differentiable and the domain of \(f\) is \(\mathbb{R}^d\) and \(f\) is \(L\)-smooth. Then the gradient descent with a fixed step size \(\gamma\leq 2/L\) satisfies \[|f(x^{(T)}) - f^*|\leq \frac{2\|x^{(0)} - x^*\|^2_2}{\gamma(2-\gamma L) }\frac{1}{T}.\] In the case of backtracking, we will replace \(\gamma\) by \(\beta/L\).
Define \(\delta_t = f(x^{(t)})-f^*\). Then subtracting \(f^*\) from both sides of equation \(\eqref{eq::lem3}\), we obtain \[\begin{align*} \delta_{t+1}&\leq \delta_t - \frac{\gamma}{2}(2-\gamma L)\|\nabla f(x^{(t)})\|^2\\ &\leq \delta_t - \frac{\gamma}{2}(2-\gamma L) \frac{|f(x^{(t)})-f^*|^2}{ \|x^{(0)}-x^*\|^2}\qquad \mbox{by \eqref{eq::lem5}}\\ &= \delta_t - \frac{\gamma}{2}(2-\gamma L) \frac{\delta_t^2}{ \|x^{(0)}-x^*\|^2}\\ &= \delta_t - \eta \delta_t^2, \end{align*}\] where \(\eta =\frac{ \gamma(2-\gamma L) }{2 \|x^{(0)}-x^*\|^2}>0\).
To obtain the rate of \(\delta_t\), we use the idea of telescoping. Dividing both sides in the above inequality by \(\frac{1}{\delta_{t+1}\delta_t}\) leads to \[\frac{1}{\delta_t} \leq \frac{1}{\delta_{t+1}} - \eta\frac{\delta_t}{\delta_{t+1}}\] and after rearrangements, this becomes \[\eta\frac{\delta_t}{\delta_{t+1}}\leq \frac{1}{\delta_{t+1}}-\frac{1}{\delta_{t}}.\] Equation \(\eqref{eq::lem4}\) shows that \(\delta_{t+1}<\delta_t\), so \(\frac{\delta_t}{\delta_{t+1}}>1\), which further implies \[\eta < \eta\frac{\delta_t}{\delta_{t+1}}\leq \frac{1}{\delta_{t+1}}-\frac{1}{\delta_{t}}.\] By summing the above inequality for \(t=0\) to \(t=T-1\), we obtain \[\eta T < \frac{1}{\delta_T} - \frac{1}{\delta_0} <\frac{1}{\delta_T}.\] Thus, we conclude \[f(x^{(T)})-f^* = \delta_T < \frac{1}{\eta T} = \frac{2\|x^{(0)}-x^*\|^2}{\gamma (2-\gamma L)} \frac{1}{T} = O(T^{-1}).\]
In the case of asymptotic analysis, we will say that gradient descent has a convergence rate \(O(1/t)\). Or we say that the gradient descent algorithm finds an \(\epsilon\)-suboptimal point in \(O(1/\epsilon)\) iterations.
If we further assume that the function not only has a Lipschitz gradient but is also \(m\)-strongly convex, i.e., \[f(x) - \frac{m}{2} \|x\|^2_2\] is convex for some \(m>0\), then the gradient descent algorithm converges even faster.
Theorem 3. Assume that \(f\) is convex and differentiable and the domain of \(f\) is \(\mathbb{R}^d\) and the gradient \(\nabla f(x)\) is Lipschitz continuous with a constant \(L>0\). Moreover, \(f\) is \(m-\)strongly convex for some \(m>0\). Suppose that \(f^* = f(x^*)= 0\). Then the gradient descent with a fixed step size \(\gamma\leq 1/L\) satisfies \[%|f(x^{(t)}) - f^*|\leq \|x^{(t)} - x^*\|^2_2\leq (1-\gamma m)^t %c^t \frac{L}{2} \|x^{(0)} - x^*\|^2_2,\quad f(x^{(t)}) \leq f(x^{(0)})\left(1-\frac{m^2}{L^2} + m^2\left(\gamma-\frac{1}{L}\right)^2\right)^{t} %\frac{L}{2}(1-\gamma m)^{t}\|x^{(0)} - x^*\|_2^2\] Note that the fastest choice of \(\gamma\) is \(\gamma=1/L\), leading to \[\|x^{(t)} - x^*\|_2\leq (1-m/L)^{t/2} \|x^{(0)} - x^*\|_2.\]
Convergence in distance. The proof relies on the following two facts.
Fact 1: When \(f(x)\) is \(m\)-strongly convex, we have \[f(x^*) - f(x^{(t)}) \geq (x^*-x^{(t)})^T \nabla f(x^{(t)}) + \frac{m}{2}\|x^*-x^{(t)}\|^2,\] which implies \[\begin{equation} (x^{(t)}-x^*)^T \nabla f(x^{(t)}) \geq f(x^{(t)})- f(x^*) + \frac{m}{2}\|x^*-x^{(t)}\|^2. \label{eq::pf1} \end{equation}\]
Fact 2: because the gradient is \(L\)-Lipschitz and \(f(x)\geq f(x^*) = 0\), for any point \(x\) we have \[\begin{align*} -f(x) &\leq f\left(x-\frac{1}{L}\nabla f(x)\right) - f(x)\\ &\overset{\eqref{eq::lem1}}{\leq} -\frac{1}{L}\|\nabla f(x)\|^2 + \frac{1}{2L}\|\nabla f(x)\|^2 \\ & = -\frac{1}{2L}\|\nabla f(x)\|^2. \end{align*}\] Thus, \[\begin{equation} f(x) \geq \frac{1}{2L}\|\nabla f(x)\|^2. \label{eq::pf2} \end{equation}\]
Consider the step at \(t+1\). Then \[\begin{align*} \|x^{(t+1)}-x^*\|^2 &=\|x^{(t)} - x^* - \gamma\nabla f(x^{(t)})\|^2 \\ &= \|x^{(t)}-x^*\|^2 - 2\gamma (x^{(t)}-x^*)^T \nabla f(x^{(t)}) + \gamma^2 \|\nabla f(x^{(t)})\|^2\\ &\overset{\eqref{eq::pf1}}{\leq } \|x^{(t)}-x^*\|^2 -2\gamma(f(x^{(t)})- f(x^*) + \frac{m}{2}\|x^*-x^{(t)}\|^2) + \gamma^2 \|\nabla f(x^{(t)})\|^2\\ & = \|x^{(t)}-x^*\|^2 (1- \gamma m) - 2\gamma f(x^{(t)}) + \gamma^2 \|\nabla f(x^{(t)})\|^2\\ %&\overset{\eqref{eq::pf2}}{\leq } \|x^{(t)}-x^*\|^2 (1- \gamma m) - 2\gamma\frac{1}{2L}\|\nabla f(x^{(t)})\|^2 + \gamma^2 \|\nabla f(x^{(t)})\|^2-2\gamma f(x^{(t)})(1-\gamma L)\\ &\overset{\eqref{eq::pf2}}{\leq } \|x^{(t)}-x^*\|^2 (1- \gamma m) - 2\gamma\frac{1}{2L}\|\nabla f(x^{(t)})\|^2 + \gamma^2 \|\nabla f(x^{(t)})\|^2\\ & = \|x^{(t)}-x^*\|^2 (1- \gamma m) - \frac{\gamma}{L} (1-\gamma L)\|\nabla f(x^{(t)})\|^2\\ &\leq \|x^{(t)}-x^*\|^2 (1- \gamma m) \end{align*}\] whenever \(\gamma<\frac{1}{L}\). Note that the quantity \(\frac{\gamma}{L} (1-\gamma L)\|\nabla f(x^{(t)})\|^2\) is the minimal improvement that the method will make per each iteration.
Convergence of objective function. Before proving the convergence of objective function, we note another fact about gradient descent.
Fact 3: The \(m\)-strongly convexity implies \(\|\nabla f(x)\|\geq m \|x-x^*\|\) and the \(L\)-Lipschitz gradient implies \(f(x)\leq \frac{L}{2} \|x-x^*\|^2\). Thus, \[\begin{equation} \|\nabla f(x)\|^2 \geq m^2 \|x-x^*\|^2 \geq \frac{2m^2}{L} f(x). \label{eq::pf3} \end{equation}\]
Using the fact that the gradient is \(L\)-Lipschitz, we have \[f(x) -f (y) \leq (x-y)^T \nabla f(y) + \frac{L}{2}\|x-y\|^2.\] Thus, \[\begin{align*} f(x^{(t+1)}) - f(x^{(t)})& =f(x^{(t)} - \gamma \nabla f(x^{(t)})) - f(x^{(t)})\\ &\leq -\gamma \|\nabla f(x^{(t)})\|^2 + \frac{L\gamma^2}{2} \|\nabla f(x^{(t)})\|^2\\ & = -\gamma\left(1-\frac{L\gamma}{2}\right)\|\nabla f(x^{(t)})\|^2\\ &\overset{\eqref{eq::pf3}}{\leq} -\gamma\left(1-\frac{L\gamma}{2}\right)\frac{2m^2}{L} f(x^{(t)}). \end{align*}\] Thus, after rearrangements, \[\begin{align*} f(x^{(t+1)}) &\leq f(x^{(t)})\left(1-2\gamma \frac{m^2}{L}\left(1-\frac{L\gamma}{2}\right)\right)\\ & = f(x^{(t)}) \left(1- \frac{2m^2}{L}\gamma + m^2\gamma^2\right)\\ & = f(x^{(t)}) \left(1-\frac{m^2}{L^2} + m^2\left(\gamma-\frac{1}{L}\right)^2\right) \end{align*}\] By telescoping, we conclude \[f(x^{(t)}) = f(x^{(0)})\left(1-\frac{m^2}{L^2} + m^2\left(\gamma-\frac{1}{L}\right)^2\right)^t,\] which completes the proof.
The key improvement of being strongly convex is the Fact 3, which provides a lower bound on the gradient through distance \(\|x-x^*\|\). This allows us to establish a quantity like \(f(x^{(t+1)})-f^* = \delta_{t+1} \leq (1-\omega) \delta_t\) (\(\omega\) is some constant), which leads to a linear (exponential) convergence. If \(f\) is only convex, then we can only bound the gradient from below using equation \(\eqref{eq::lem5}\), which leads to \(\delta_{t+1} \leq \delta_t (1-\eta \delta_t)\), so we can only establish a polynomial convergence.
With the strongly convex assumption, the convergence rate is \(O((1-m/L)^{t})\), exponentially fast! Thus, finding an \(\epsilon\)-suboptimal point only takes \(O(\log(1/\epsilon))\) iterations. The quantity \(L/m\) is called the condition number. A higher condition number implies a slower convergence rate.
When the convergence rate is an exponential function with respect to \(t\), we say that this algorithm has a linear convergence rate. This is because it looks like a linear on a semi-log plot (log error versus iterations). The following is the log error versus iterations in the previous example:
Finally, we conclude this by introducing an interesting result about the limitation of convergence rate. The first-order method is the collection of iterative methods for updating \(x^{(t)}\) using the value of \(f(x)\) and its gradients \(\nabla f(x)\). The gradient descent is a first-order method.
Theorem 4 (Nesterov). For any \(t\leq (d-1)/2\) and any starting point \(x^{(0)}\), there exists \(f\) , a convex, differentiable function with Lipschitz gradient, such that any first-order method satisfies \[f(x^{(t)}) - f^* \geq \frac{3L \|x^{(0)}-x^*\|^2_2}{32 (t+1)^2}.\]
Namely, the convergence rate cannot be faster than \(O(1/t^2)\). Note that the gradient descent approach achieve the exponential rate due to the strongly convex assumption. This lower bound can be achieved by using the Nesterov accelerated gradient descent algorithm10, a modification from the gradient descent algorithm.
Remark.
Similar to the EM algorithm, the gradient descent also suffers from local minima. Note that when we assume \(f\) to be a convex function, there is only one unique minimum. If \(f\) is non-convex and has multiple local minima, one has to choose the initial point multiple times and apply the gradient descent algorithm with different starting points to increase the chance of obtaining the global minimum.
If you are interested in the proof or related topics, I would recommend reading the followings:
https://perso.telecom-paristech.fr/rgower/pdf/M2_statistique_optimisation/grad_conv.pdf
Introductory Lectures on Convex Programming, Volume I: Basic course by Yu. Nesterov.
In analyzing property of gradient descent, often we can consider the continuous limit of the gradient descent. When the step size \(\gamma\approx 0\), the gradient descent behaves like the a gradient flow \(\eta_x(t)\) with \[\eta_x(0) = x, \quad \frac{\partial}{\partial t}\eta_x(t) = -\nabla f(\eta_x(t))\] and \(x = x^{(0)}\) is the initial point. Thus, analyzing the property of such a gradient flow gives us insight about the property of gradient descent. This idea is particularly useful when the function \(f\) is non-convex and may contain multiple critical points11. A very important assumption to obtain good gradient flow is that \(f\) is a Morse function, meaning that all critical points are well-separated.
The analysis on the continuous-time gradient flow gives us some insight on the gradient descent/ascent algorithm. The study of the continuous-time gradient flow is a common topic in dynamical systems, a branch in applied mathematics. For a gradient system, the collection \(\mathcal{C}= \{x: \nabla f(x) = 0\}\) is called stationary points or critical points. As can be seen easily, the gradient flow will stop when it reaches any point within \(\mathcal{C}\). So for any element of \(\mathcal{C}\), we can define its basin of attraction: \[D(c) = \{x: \eta_x(\infty)\}\] Interestingly, this implies that the gradient decent may not reach a local maximum–it may end up being at a saddle point or a local minimum. Here are two interesting properties characterizing the gradient flow and stationary points.
Theorem 5. Let \(\eta_x{(\infty)} = \lim_{t\rightarrow \infty}\eta_x{(t)}\in \mathcal{C}\) be the destination of a gradient descent flow starting at \(x\). Let \((\lambda^+_{\min}(x), v^+_{\min}(x))\) be the smallest positive eigenvalue and its eigenvector of \(\nabla \nabla f(x)\).
Terminal flow: If \(\lambda^+_{\min}(\eta_x{(\infty)})\) is unique, then the terminal direction of the gradient flow \[\lim_{t\rightarrow\infty} \frac{\eta_x{(\infty)} - \eta_x{(t)}}{\|\eta_x{(\infty)} - \eta_x{(t)}\|} = - v^+_{\min}(\eta_x{(\infty)}).\]
Stable manifold theorem: If \(f\) is a Morse function (all critical points are well-separated), the basin of attraction \(D(c)\) forms a \(k\)-dimensional manifold if there are \(k\) positive eigenvalues of \(\nabla\nabla f(c)\).
Recall that in the likelihood inference, we are often interested in finding the MLE \[\hat\theta = {\sf argmax}_{\theta} \ell(\theta|X_1,\cdots, X_n),\] where \(\ell(\theta|X_1,\cdots, X_n) = \frac{1}{n}\sum_{i=1}^n \ell(\theta|X_i)\) is the log-likelihood function.
As you can see, the log-likelihood function behaves like a sample mean. In fact, it can be viewed as an estimator of the population likelihood function \[\ell(\theta) = \mathbb{E}(\ell(\theta|X_i))\] and the MLE is estimating the population MLE \[\theta^* = {\sf argmax}_\theta \ell(\theta).\]
In establishing the asymptotic normality of the MLE, we often assume that the second derivative (known as the Hessian matrix) \(H(\theta) = \nabla \nabla \ell(\theta)\) exists and is smooth so that we can apply a Taylor expansion to show that \[\sqrt{n}(\hat\theta - \theta^*)\approx N(0, \Sigma),\] where \[\Sigma = H^{-1}(\theta^*) I(\theta^*)H^{-1}(\theta^*),\quad I(\theta) = \mathbb{E}(\nabla \ell(\theta) \nabla \ell(\theta)^T).\] We can estimate \(\Sigma\) via a plug-in estimate \[\hat\Sigma = \hat H^{-1}(\hat\theta) \hat I(\hat\theta)\hat H^{-1}(\hat\theta),\quad\hat H(\theta) = \frac{1}{n}\sum_{i=1}^n \nabla\nabla \ell(\theta|X_i),\quad \hat I(\theta) = \frac{1}{n}\sum_{i=1}^n \nabla \ell(\theta|X_i)\nabla \ell(\theta|X_i)^T\] and use it to construct a confidence interval. The estimator \(\hat\Sigma\) is known as the sandwich estimator and it is a consistent estimator of \(\Sigma\) if \(H(\theta)\) is smooth. Note that you may find \(\Sigma\) somewhat different from the conventional form. This is because the current form of \(\Sigma\) does not require the parametric model to be correct. If the parametric model is correct, \(\Sigma\) will reduce to \(I^{-1}(\theta^*)\), the usual Fisher information matrix.
Here is an interesting relation from the asymptotic normality to the optimization problem we just described. To establish the asymptotic normality, clearly we need the matrix \(H(\theta^*)\) to be invertible and it has to be smooth. The invertibility of a matrix is related to its eigenvalues– a matrix is invertible if and only if its eigenvalues are non-zero. In particular, since \(\theta^*\) is the maximizer, we need all eigenvalues of \(H(\theta^*)\) to be negative. This is essentially saying that we need \(H(\theta)\) to be be strongly concave at the neighborhood of \(\theta^*\) to establish the asymptotic normality and validity of a normal confidence interval.
To sum up, we have the following beautiful theorems connecting statistical validity and optimization feasibility.
Theorem 6. Suppose that the (population) likelihood function \(\ell(\theta)\) satisfies the following conditions:
\(\theta^*\) is the unique maximizer of \(\ell(\theta)\) in the interior of the parameter space \(\Theta\).
\(H(\theta) = \nabla \nabla \ell(\theta)\) is smooth.
all eigenvalues of \(H(\theta^*)\) are negative.
Suppose that we apply a gradient ascent algorithm to find the MLE: \[\hat\theta^{(t+1)} =\hat\theta^{(t)} +\gamma \nabla \ell(\hat\theta^{(t)}|X_1,\cdots, X_n).\]
Then we have the following statistical validity:
Consistency of the MLE: \(\hat \theta \overset{P}{\rightarrow} \theta^*.\)
Asymptotic normality of the MLE: \(\sqrt{n}(\hat\theta - \theta^*)\overset{d}{\rightarrow} N(0, \Sigma).\)
Validity of sandwich estimator: \(\hat \Sigma \overset{P}{\rightarrow} \Sigma.\)
as well as the following optimization feasibility: there exists \(r_0>0\) such that
Validity of the gradient ascent: if \(\hat\theta^{(0)}\in B(\theta^*,r_0)\), then \(\hat\theta^{(\infty)} = \lim_{t\rightarrow\infty}\hat\theta^{(t)} = \hat\theta\) when the step size \(\gamma\) is sufficiently small.
Linear convergence of gradient ascent: \(\|\hat\theta^{(t)} - \hat\theta\|\leq (1-\omega)^{t/2} \|\hat\theta^{(0)} - \hat\theta\|\) for some \(\omega\in (0,1)\).
The beauty of the above theorem is that it provides a direct association between statistical stability and algorithmic convergence. Note that in the above theorem, we only require local strong concavity of \(\ell(\theta)\) around \(\theta^*\). If the entire likelihood function is strongly concave, then the requirement \(\theta^{(0)} \in B(\theta^*,r_0)\) can be relaxed to \(\theta^{(0)} \in \Theta\)–any initial point will work. In fact, a similar theorem can be applied to other statistical estimators as long as they can be written in the form of ‘maximizing some sample mean objective function’. A class of this type of statistical problem is called an M-estimator.
In a sense, Theorem 6 establishes the computability of Frequentist inference: we can find the MLE using the gradient ascent method under proper conditions. This, together with what we have learned in the rejection sampling that we draw a proposal from the prior distribution and reject it with the probability \(P(\mbox{Reject}) = \frac{L(\theta|X_1,\cdots,X_n)}{L(\hat\theta|X_1,\cdots, X_n)}\), implies computability of a Bayesian inference: we can draw from the posterior distribution by the rejection sampling with the computed MLE. This allows us to handle problems where there is no closed-form MLE (one famous example is the logistic regression).
The stochastic gradient descent (SGD) is a popular approach to perform minimization of a function \(f(x) = \frac{1}{n}\sum_{i=1}^n f_i(x)\). There are many concrete example of such case; for instance, to find the MLE, we are trying to minimize the negative log-likelihood function (here we define \(\ell_n\) as the log-likelihood function divided it by \(n\)): \[-\ell_n(\theta) = -\frac{1}{n}\log \prod_{i=1}^n p(X_i;\theta) = -\frac{1}{n}\sum_{i=1}^n \log p(X_i;\theta) =- \frac{1}{n}\sum_{i=1}^n \ell(\theta|X_i).\] The traditional gradient descent will start with an initial point \(\theta^{(0)}\) and update it according to \[\theta^{(t+1)} \leftarrow \theta^{(t)} + \gamma_t \nabla \ell_n(\theta^{(t)}).\] Note that here we allow the step size \(\gamma=\gamma_t\) to be dependent on the iteration.
Let \(x^{(0)}\) be the initial point and \(\gamma_t\) be a sequence of step size. To minimizes the function \(f(x) = \frac{1}{n}\sum_{i=1}^n f_i(x)\), the stochastic gradient descent (SGD) uses the following update \[x^{(t+1)} \leftarrow x^{(t)} + \gamma_t \nabla f_{i_t} (x^{(t)}),\] where \(i_t \in \{1,2,\cdots, n\}\) is a chosen index at iteration \(t\). Two common rules for choosing \(i_t\) are
Randomized rule: \(i_t\) is a uniform distribution over \(\{1,2,\cdots, n\}\).
Cyclic rule: choose \(i_t\) as \(1,2,3,\cdots, m, 1,2,3,\cdots\).
Randomized rule is more popular in practice. The idea of randomized rule is that \[\mathbb{E}(\nabla f_{i_k}(x)) = \nabla f(x),\] i.e., each SGD update is an unbiased estimate of the gradient.
The trajectory \(\{x^{(0)}, x^{(1)},\cdots, \}\) forms a Markov chain. Moreover, if we fix the step size \(\gamma_t=\gamma>0\), this Markov chain is homogeneous.
The main appeals of SGD are (i) the iteration cost is independent of \(n\), the number of functions (in statistics, it is often the sample size), and (ii) it may provide big savings in terms of memory usage (no need to read the entire dataset). Thus, SGD is powerful when we have a large-scale optimization problem. Specifically, assume that computing the function \(f_i(x)\) requires one unit of computational cost. Then, one update of traditional gradient descent algorithm requires \(O(n)\) cost whereas one iteration of SGD only requires \(O(1)\) cost.
If we keep the step size \(\gamma\) the same across iterations, the SGD will not converge (it forms a Markov chain so at the end it will behave like a random point from the stationary distribution). Thus, to ensure SGD converges, a standard approach is to choose diminishing step sizes such as \(\gamma_t = 1/t\).
Convergence rate. Recall that when \(f\) has a Lipschitz gradient, traditional gradient descent has a convergence rate \(O(1/t)\). In the case of SGD, however, the convergence rate is \(O(1/\sqrt{t})\)12. More formally, \[\mathbb{E}(f(x^{(t)})) - f^* = O(1/\sqrt{t}).\] Thus, SGD has a slower convergence rate. Moreover, when \(f\) is strongly convex with a Lipschitz gradient, the gradient descent is of the rate \(O(c^t)\) for some \(c<1\) but under the same condition, SGD gives a rate \(O(1/t)\). Although it is cheap to evaluate one iteration of SGD, the convergence rate is a lot slower than the traditional gradient descent.
Does this imply that SGD is not useful? Not at all! Consider the case where the function \(f\) has a Lipschitz gradient and our goal is to find a minimizer that is at least \(10^{-3}\) close to the optimal. Assume that \(f\) is the average of \(n=10^6\) function (loosely speaking, the sample size is \(10^6\)). How much cost will it take for gradient descent and SGD to achieve this accuracy? For gradient descent, the convergence rate is \(O(1/t)\), so \(O(1/t) = 10^{-3}\Rightarrow t= O(10^3)\), which implies that we need \(O(10^3)\) iterations. Each iteration costs \(O(n) = O(10^6)\) computational budget so the total cost is \(O(nt) = O(10^9)\). For SGD, the convergence rate is \(O(1/\sqrt{t})\). Thus, we need \(t = O(10^{6})\) iterations. However, each iteration only cost \(O(1)\) budget so the total cost is \(O(t) = O(10^6)\), a \(1000\)-times saving from traditional gradient descent!
This property has made SGD very popular in large-scale computation especially in the modern big data era. A sample of size \(n=10^6\) is considered as big data set nowadays but as you can see, the larger the dataset, the more benefit SGD has over traditional gradient descent.
Mini-batch approach. A popular variant of SGD is the mini-batch approach. Instead of randomly choose one function, we randomly choose \(b\) functions and use the gradient of the average as the update. In more details, let \(I_t \subset\{1,2,\cdots, n\}\) be a set of size \(|I_t| = b<<m\). We then use the following update: \[x^{(t+1)} \leftarrow x^{(t)} - \gamma_t \frac{1}{b} \sum_{i \in I_t} f_i (x^{(t)}).\]
Averaged SGD. A new approach in SGD is the averaged SGD approach. The idea is very simple. After running the SGD and observing the trajectory \(\{x^{0},x^{1},\cdots, x^{T}\}\), instead of using the last point as the output minimizer, we use the average of the last few points as the minimizer. Namely, we choose a number \(\tau\) and then use \[\tilde{x}_{T} = \frac{1}{\tau}\sum_{k=0}^{\tau-1} x^{(T-k)}\] as the minimizer. This idea has been shown to enjoy both a fast convergence rate and a low computational cost13.
There is another set of SGD approach that uses the idea of additive noise. This is called the perturbed gradient descent14. Let \(f\) be the function that we wish to minimize. Starting with an initial point \(x^{(0)}\) and a sequence of step sizes \(\gamma_t\), we use the following update: \[x^{(t+1)}\leftarrow x^{(t)} - \gamma_t \nabla f(x^{(t)}) +\gamma_t Z_t,\] where \(Z_1,Z_2,\cdots,\) are IID mean \(0\) random noises (often we choose them to be multivariate normal). This approach does not save computational cost (because we still evaluate a full gradient) but it may escape from (shallow) local minima due to the additive noise \(Z_t\).
Note that a similar idea to the additive noise is the stochastic gradient Langevin dynamics (SGLD)15, which uses \[x^{(t+1)}\leftarrow x^{(t)} - \gamma_t \nabla f(x^{(t)}) +\sqrt{\gamma_t} Z_t,\] where \(Z_t\) is a Gaussian vector with identity covariance matrix. When the step size \(\gamma_t= \gamma\) is fixed, SGLD forms a Markov chain with a stationary distribution \(\propto e^{-\gamma f(x)}\). Thus, SGLD is often used as an alternative to MCMC although it may be used as a method to find the minimizer. Note that one can choose \(f(x) = \log \pi(\theta|X_1,\cdots, X_n)\) then SGLD allows us to sample from the posterior distribution16.
There are also other paradigms such as fiducial inference (https://en.wikipedia.org/wiki/Fiducial_inference), but the Frequentist and Bayesian paradigms are the two major ones.↩︎
see, e.g., Chapter 2 of All of nonparametric statistics by Larry Wasserman and Chapter 3.5 of The bootstrap and Edgeworth expansion by Peter Hall.↩︎
https://www.amstat.org/asa/files/pdfs/P-ValueStatement.pdf↩︎
see https://en.wikipedia.org/wiki/Jeffreys_prior for more details.↩︎
https://en.wikipedia.org/wiki/Bernstein%E2%80%93von_Mises_theorem↩︎
see, e.g., the famous paper Ghosal et al. (2000) “Convergence rates of posterior distributions": https://projecteuclid.org/euclid.aos/1016218228↩︎
See http://mitliagkas.github.io/ift6085/ift-6085-lecture-6-notes.pdf↩︎
see, e.g., https://arxiv.org/pdf/1602.04915.pdf and https://arxiv.org/abs/1807.04431.↩︎
see, .e.g., Nemirosvki et al. (2009) Robust stochastic optimization approach to stochastic programming.↩︎
See, e.g., https://www.di.ens.fr/~fbach/colt_2018_pillaud_vivien.pdf and http://auai.org/uai2018/proceedings/papers/71.pdf.↩︎
See https://arxiv.org/abs/1703.00887v1 for more details.↩︎
For applications in Statistics, see http://www.jmlr.org/papers/volume17/teh16a/teh16a.pdf↩︎
see https://www.ics.uci.edu/~welling/publications/papers/stoclangevin_v6.pdf for more details.↩︎