Bias-variance tradeoff

According to wiki:

In statistics and machine learning, the bias–variance tradeoff is the property of a set of predictive models whereby models with a lower bias in parameter estimation have a higher variance of the parameter estimates across samples, and vice versa.

Models with low bias are usually more complex (e.g. higher-order regression polynomials), enabling them to represent the training set more accurately. In the process, however, they may also represent a large noise component in the training set, making their predictions less accurate - despite their added complexity. In contrast, models with higher bias tend to be relatively simple (low-order or even linear regression polynomials) but may produce lower variance predictions when applied beyond the training set.

Bias–variance decomposition of mean squared error (MSE):

We assume $y=f\left(x\right)+\epsilon$, where $E\left(\epsilon \right)=0$ and $\text{Var}\left(\epsilon \right)={\sigma }^2$. Our goal is to find a function $\hat{f}\left(x\right)$ that makes MSE of $\hat{f}$, $E\left\{\right(y-\hat{f}{)}^T(y-\hat{f})\}$, minimum.

The Bias-Variance decomposition of MSE proceeds as follows: $$\begin{array}{ccc}& & E\left\{\right(y-\hat{f}{)}^T(y-\hat{f})\}=\{E(y-\hat{f}){\}}^{T}E(y-\hat{f})+\text{Var}(y-\hat{f})\\ & =& ||\text{Bias}\left(\hat{f}\right)|{|}^2+\text{Var}\left(y\right)+\text{Var}\left(\hat{f}\right)-2\text{cov}(y,\hat{f})\\ & =& ||\text{Bias}\left(\hat{f}\right)|{|}^2+\text{Var}\left(\hat{f}\right)+{\sigma }^2,\end{array}$$ where $$\begin{array}{ccc}\text{cov}(y,\hat{f})& =& E\left(y\hat{f}\right)-E\left(y\right)E\left(\hat{f}\right)\\ & =& E\left[\right(f+\epsilon \left)\hat{f}\right]-E(f+\epsilon )E\left(\hat{f}\right)\\ & =& fE\left(\hat{f}\right)+E\left(\epsilon \hat{f}\right)-fE\left(\hat{f}\right)\\ & =& E\left(\epsilon \hat{f}\right)\\ & =& 0,\end{array}$$ since $\epsilon \perp \hat{f}$ or they are independent. (Question : why independent implies $\epsilon \perp \hat{f}$, which implies $E\left(\epsilon \hat{f}\right)=0$).

Bias-variance tradeoff

• models including many covariates leads to have low bias but high variance.
• models including few covariates leads to high bias but low variance.

Hence, we need criteria that both take in account model complexity (number of predictors) and quality of fit.

$R^2$

Denifition: $$\begin{array}{c}{\text{R}}^2=1-\frac{RSS}{TSS}=1-\frac{{\sum }_i(y_i-\widehat{f_i}{)}^2}{{\sum }_i(y_i-\overline{y}{)}^2},\end{array}$$ where TSS is total sum of squares, RSS is residual sum of squares. And define $\text{ESS}={\sum }_i(\hat{f}-\overline{y}{)}^2$ is explained sum of squares, also called the regression sum of square. $R^2$ is based on the assumption that $TSS=RSS+ESS$. Under the linear regression model setting satisfies this assumption usually.

Proof:

$$\begin{array}{ccc}& & \sum _i(y_i-\overline{y}{)}^2=\sum (y_i-\widehat{f_i}+\widehat{f_i}-\overline{y}{)}^2\\ & =& \sum _i(y_i-\widehat{f_i}{)}^2+\sum _i(\widehat{f_i}-\overline{y}{)}^2+2\sum _i(y_i-\widehat{f_i})(\widehat{f_i}-\overline{y})\\ & =& RSS+ESS+2\sum _i{\hat{e}}_i(\widehat{f_i}-\overline{y})(\widehat{f_i}=\widehat{y_i}=X\hat{\beta }\text{in linear model})\\ & =& RSS+ESS+2\sum _i{\hat{e}}_i(\widehat{y_i}-\overline{y})\\ & =& RSS+ESS+2\sum _i{\hat{e}}_i\widehat{y_i}-2\overline{y}\sum _i{\hat{e}}_i\end{array}$$ Then, the reamining part is to prove ${\sum }_i{\hat{e}}_i(\widehat{y_i}-\overline{y})=0$.

Firstly, ${\sum }_i{\hat{e}}_i\widehat{y_i}=e^{T}Hy=y^T(I-H)Hy=0$ due to $H$ idempotent. Then if we can show ${\sum }_i{\hat{e}}_i=0$, our proof is done. However, this can not be shown for a model without an intercept.

Ajusted $R^2$¹

Define adusted $R^2$ as

$R_a^2=1-\frac{RSS/d{f}_e}{TSS/d{f}_t}=1-\frac{RSS/(n-p-1)}{TSS/(n-1)}$ = $1-\frac{Var\left({\hat{\sigma }}^2\right)}{Va{r}_{tot}}$, where $Var\left({\hat{\sigma }}^2\right)$ is the sample variance of the estimated residuals and $Va{r}_{tot}$ is sample variance of dependent variable. They can be considered as unbiased estimates of the population variances of the errors and of the dependent variable. Hence, adjusted $R^2$ can be interpreted as an unbiased estimator of the population $R^2$.

Laplace’s approximation

Define an index set of the active predictors $\gamma =\{j:{\beta }_j\ne 0\}$ for $j=1,2,\dots ,p$. So $\gamma$ can be treated as a model we consider. The Bayesian approach to the model selection is to maximize the posterior distribution of a model given the data $\pi \left(\gamma |y\right)$. By Bayes theorem, $\pi \left(\gamma |y\right)\propto f\left(y|\gamma \right)\pi \left(\gamma \right)$, where $$f\left(y|\gamma \right)=\int f(y|{\theta }_{\gamma },\gamma )\pi \left({\theta }_{\gamma }|\gamma \right)d{\theta }_{\gamma }=\int \exp \left[\log \left(f\right(y|{\theta }_{\gamma },\gamma \left)\pi \right({\theta }_{\gamma }|\gamma \left)\right)\right]d{\theta }_{\gamma }.$$ We expand $\log \left(f\right(y|{\theta }_{\gamma },\gamma \left)\pi \right({\theta }_{\gamma }|\gamma \left)\right)$ by Taylor expansion about its posterior mode ${\hat{\theta }}_{\gamma }$ where $f(y|{\theta }_{\gamma },\gamma )\pi \left({\theta }_{\gamma }|\gamma \right)$ attains the maximimum. Thus, $$\log \left(f\right(y|{\theta }_{\gamma },\gamma \left)\pi \right({\theta }_{\gamma }|\gamma \left)\right)\approx \log \left(f\right(y|{\hat{\theta }}_{\gamma },\gamma \left)\pi \right({\hat{\theta }}_{\gamma }|\gamma \left)\right)+({\theta }_{\gamma }-{\hat{\theta }}_{\gamma }){\nabla }_{\theta }Q{|}_{\hat{\theta }}+\frac{1}{2}({\theta }_{\gamma }-{\hat{\theta }}_{\gamma }{)}^T{H}_{{\hat{\theta }}_{\gamma }}({\theta }_{\gamma }-{\hat{\theta }}_{\gamma }),$$ where $Q_{{\hat{\theta }}_{\gamma }}=\log \left(f\right(y|{\theta }_{\gamma },\gamma \left)\pi \right({\theta }_{\gamma }|\gamma \left)\right)$ and $H_{{\hat{\theta }}_{\gamma }}$ is a Hessian matrix such that $H_{{\hat{\theta }}_{\gamma }}=\frac{{\partial }^{2}Q}{\partial {\theta }^T\partial \theta }{|}_{{\hat{\theta }}_{\gamma }}$. Note that ${\nabla }_{\theta }{Q}_{\hat{\theta }}=0$ and $H_{{\hat{\theta }}_{\gamma }}$ is negative definite at ${\hat{\theta }}_{\gamma }$.

Therefore, $$\begin{array}{ccc}f\left(y|\gamma \right)& =& \int \exp \left[\log \left(f\right(y|{\theta }_{\gamma },\gamma \left)\pi \right({\theta }_{\gamma }|\gamma \left)\right)\right]d{\theta }_{\gamma }\\ & \approx & \int \exp \{Q_{{\hat{\theta }}_{\gamma }}-\frac{1}{2}({\theta }_{\gamma }-{\hat{\theta }}_{\gamma }{)}^T{\tilde{H}}_{{\hat{\theta }}_{\gamma }}({\theta }_{\gamma }-{\hat{\theta }}_{\gamma })\}d{\theta }_{\gamma }\\ & =& \exp \left(Q_{{\hat{\theta }}_{\gamma }}\right)\int \exp \{-\frac{1}{2}({\theta }_{\gamma }-{\hat{\theta }}_{\gamma }{)}^T{\tilde{H}}_{{\hat{\theta }}_{\gamma }}({\theta }_{\gamma }-{\hat{\theta }}_{\gamma })\}d{\theta }_{\gamma }\\ & =& f(y|{\hat{\theta }}_{\gamma },\gamma )\pi \left({\hat{\theta }}_{\gamma }|\gamma \right){\left|2\pi {\tilde{H}}_{{\hat{\theta }}_{\gamma }}^{-1}\right|}^{\frac{1}{2}}\\ & =& f(y|{\hat{\theta }}_{\gamma },\gamma )\pi \left({\hat{\theta }}_{\gamma }|\gamma \right)\frac{(2\pi {)}^{|\gamma |/2}}{|{\tilde{H}}_{{\hat{\theta }}_{\gamma }}{|}^{1/2}},\end{array}$$ where ${\tilde{H}}_{{\hat{\theta }}_{\gamma }}=-H_{{\hat{\theta }}_{\gamma }}$ is the negative and Hessian matrix and the observed Fisher information matrix. Taking log of it, we obtain

$$\begin{array}{}2\log f\left(y|\gamma \right)=2\log f(y|{\hat{\theta }}_{\gamma },\gamma )+2\log \pi \left({\hat{\theta }}_{\gamma }|\gamma \right)+|\gamma |\log \left(2\pi \right)+\log |{\tilde{H}}_{{\hat{\theta }}_{\gamma }}{|}^{-1}& \text{(1)}\end{array}$$

Flat Prior and the Weak Law of Large Numbers

If we set a flat prior on the ${\theta }_{\gamma }$ such that $\pi \left({\theta }_{\gamma }\right)=1$, then each element in the matrix ${\tilde{H}}_{{\hat{\theta }}_{\gamma }}$ is $$\begin{array}{ccc}{\tilde{H}}_{mn}& =& -\frac{{\partial }^2\log f(y|{\theta }_{\gamma },\gamma )}{\partial {\theta }_m\partial {\theta }_n}{|}_{\theta ={\hat{\theta }}_{\gamma }}\\ & =& -\frac{{\partial }^2\log \left({\prod }_{i=1}^{n}f\right(y_i|{\theta }_{\gamma },\gamma \left)\right)}{\partial {\theta }_m\partial {\theta }_n}{|}_{\theta ={\hat{\theta }}_{\gamma }}\\ & =& -\frac{{\partial }^2\left({\sum }_{i=1}^n\log f\right(y_i|{\theta }_{\gamma },\gamma \left)\right)}{\partial {\theta }_m\partial {\theta }_n}{|}_{\theta ={\hat{\theta }}_{\gamma }}\\ & =& -\frac{{\partial }^2\left(\frac{1}{n}{\sum }_{i=1}^{n}n\log f\right(y_i|{\theta }_{\gamma },\gamma \left)\right)}{\partial {\theta }_m\partial {\theta }_n}{|}_{\theta ={\hat{\theta }}_{\gamma }}\end{array}$$

Since $y_1,y_2,\dots ,y_n$ are i.i.d, according to the weak law of large numbers on the random variable $z_i=n\log f(y_i|{\theta }_{\gamma },\gamma )$, we get

$$\frac{1}{n}\sum _{i=1}^{n}n\log f(y_i|{\theta }_{\gamma },\gamma )\stackrel{p}{\to }E\left[n\log f\right(y_i|{\theta }_{\gamma },\gamma \left)\right]$$ Therefore, each entry in the observed Fisher information matrix is

$$\begin{array}{ccc}{\tilde{H}}_{mn}& =& -\frac{{\partial }^2\left(E\right[n\log f(y_i|{\theta }_{\gamma },\gamma )\left]\right)}{\partial {\theta }_m\partial {\theta }_n}{|}_{\theta ={\hat{\theta }}_{\gamma }}\\ & =& -n\frac{{\partial }^2\left(E\right[\log f(y_i|{\theta }_{\gamma },\gamma )\left]\right)}{\partial {\theta }_m\partial {\theta }_n}{|}_{\theta ={\hat{\theta }}_{\gamma }}\\ & =& n{I}_{mn}\end{array}$$ for $i=1,2,\dots ,n$. So ${\tilde{H}}_{{\hat{\theta }}_{\gamma }}=n{I}_{{\hat{\theta }}_{\gamma }}$, where $I_{{\hat{\theta }}_{\gamma }}=\frac{{\partial }^{2}E\left[\log f\right(y_i|{\theta }_{\gamma },\gamma \left)\right]}{\partial {\theta }^T\partial \theta }{|}_{\theta =\hat{\theta }}$ is the Fisher information matrix for a single data point $y_i$. Thus, $|{\tilde{H}}_{{\hat{\theta }}_{\gamma }}|=n^{|\gamma |}|I_{{\hat{\theta }}_{\gamma }}|$. We plug it to the (1) and get $$\begin{array}{c}2\log f\left(y|\gamma \right)=2\log f(y|{\hat{\theta }}_{\gamma },\gamma )+2\log \pi \left({\hat{\theta }}_{\gamma }|\gamma \right)+|\gamma |\log \left(2\pi \right)-|\gamma |\log n-\log |I_{{\hat{\theta }}_{\gamma }}|.\end{array}$$ For large $n$, the item without $n$ can be negleted, and we obtain a simplified form after taking a negative sign: $$\begin{array}{}-2\log f\left(y|\gamma \right)\approx -2\log f(y|{\hat{\theta }}_{\gamma },\gamma )+|\gamma |\log n.& \text{(2)}\end{array}$$ The right-hand side of (2) is the BIC estimate for the model $\gamma$.

Laplace’s method

$${\int }_a^b{e}^{Mf\left(x\right)}\approx \sqrt{\frac{2\pi }{M|f^{{}^{\prime \prime }}\left(x_0\right)|}}{e}^{Mf\left(x_0\right)}\text{as}M\to \infty ,$$ where $f\left(x\right)$ is twice-differentiable function with a global maximum at $x_0$, $M$ is a large number, and the endpoints $a$ and $b$ could possibly be infinite.