# AIC and likelihood function applied to a dataset

First Part

I am trying to have a better understanding of what is a likelihood function and what is AIC.

Here is my current understanding: If we toss a coin 11 times and get 5 heads and 6 tails, the likelihood function is $$L(p) = P(data|p) = p^5(1-p)^6$$ p being the hyopthesis of the function which is the probability to get head.

the Maximum Likelihood is the value of p that satisfies this equation: $$\frac{dL(p)}{dp} = 5p^4(1-p)^6-6p^5(1-p)^5=0$$ Therefore, $\hat p=\frac{5}{11}$

The AIC is given by: $$AIC=2k-2ln(L)$$ k being the number of parameters (one in our example). Therefore: $$AIC=2*1-2*ln(0^5*(1-0)^6) = 2-2*ln(L)=2$$ I guess I missunderstand something! Is there a mistake? Where is it?

Second Part

I am also wondering how does this AIC calculations apply when performing regressions. For example we have this simple set of data (writing in R format):

respond_variable = c(1,2.2,3.1,4)
predictor_variable = c(1,2,6,10)


We perform a linear and a quadratic model on them. How is the AIC calculated?

Thanks a lot for your help!

Your AIC calculation has left out the dependence on $p$. It should be

$$\begin{split} AIC(p) &= 2 \times 1 - 2 \times \log \{p^5 \times (1-p)^6 \}\\ &= 2 \times 1 - 2 \times \log \{ (5/11)^5 \times (6/11)^6 \} \\ &= 2 - 2 \times 5 \log 5/11 - 2 \times 6 \log 6/11 \\ &= 17.1582 \end{split}$$

It's also left out the binomial coefficient, but this doesn't matter since it's constant for the data, and so won't affect the choice of $p$.

The AIC works the same way for regression models, only now it's a function of the vector of coefficients $\beta$. So you maximise your log-likelihood given how ever many coefficients you have, then subtract (twice) the length of $\beta$. This adjusts for the greater complexity of models with more regression parameters.

ETA: Assume we're dealing with ordinary linear regression, so we can write the model as

$$Y_i = \beta_0 + \sum_{j=1}^p \beta_jx_{ij} + \mathcal{N}(0, \sigma^2).$$

Notice that this is the same as saying that each of the $Y$'s is a normally distributed random variable whose mean is a function of the predictors $x_1 \dots x_p$ and coefficients $\beta_1 \dots \beta_p$, and whose variance is $\sigma^2$. We can thus write the likelihood as

$$L(\beta, \sigma^2) = \prod_{i=1}^n \phi\left(Y_i; \, \beta_0 + \sum_{j=1}^p \beta_jx_{ij}, \, \sigma^2 \right)$$

where $\phi$ is the density function of the normal distribution. We can also write this in a more compact way, using matrix-vector notation, as

$$\prod_i \phi(Y_i; \beta'\mathbb{x}_i, \sigma^2)$$

The log-likelihood is thus

$$\begin{split} \mathcal{L}(\beta, \sigma^2) &= \sum_i \log \phi(Y_i; \beta'\mathbb{x}_i, \sigma^2) \\ &= \sum_i \log \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left\{- \frac{(Y_i - \beta'\mathbb{x}_i)^2}{2 \sigma^2} \right\} \right) \\ &\sim -\frac{1}{2\sigma^2} \sum_i (Y_i - \beta'\mathbb{x}_i)^2 - \frac{n}{2}\log\sigma^2 \end{split}$$

and the AIC when this is maximised as a function of $\beta$ (and $\sigma^2$) is

$$AIC(\hat \beta, \hat \sigma^2) = -2 \mathcal{L}(\hat \beta, \hat \sigma^2) + 2 (p + 2)$$

The $+2$ is because we are also estimating the residual variance $\sigma^2$ and intercept term $\beta_0$. The number of fitted parameters is therefore two more than the number of predictor variables. But this can generally be ignored, since you're (almost) always estimating the residual variance and intercept from the data.

This can actually be simplified further, because the maximum likelihood estimate $\hat \sigma^2$ is just $\sum_i (Y_i - \hat\beta'\mathbb{x}_i)^2/n$ (note the denominator is NOT $n - p$), so we have

$$AIC \sim n \log \hat\sigma^2 + 2p = n \log \sum \frac{(Y_i - \hat\beta'\mathbb{x}_i)^2}{n} + 2p$$

The numerical values you get from the above formula will be different to those from the expression involving the normal density, since a bunch of arbitrary constants have been left out. That's okay, since it's the difference in AIC between models that counts, and these differences will be unaffected.

In R:

# the data -- did you interchange x and y by mistake?
y = c(1,2.2,3.1,4)
x = c(1,2,6,10)

# modelling with linear and quadratic terms
m1 <- lm(y ~ x)
m2 <- lm(y ~ x + I(x^2))

# fitted values
f1 <- fitted(m1)
f2 <- fitted(m2)

# maximum likelihood estimates for sigma -- the denominator is n, NOT n-p
s1 <- sqrt(sum((y - f1)^2)/4)
s2 <- sqrt(sum((y - f2)^2)/4)

# getting the log-likelihoods and AIC manually
L1 <- sum(dnorm(y, f1, s1, log=TRUE))
L2 <- sum(dnorm(y, f2, s2, log=TRUE))
A1 <- -2 * L1 + 2 * 3
A2 <- -2 * L2 + 2 * 4

# compare with built-in functions
l1 <- logLik(m1)
l2 <- logLik(m2)
a1 <- AIC(m1)
a2 <- AIC(m2)

a1
[1] 8.608357
A1
[1] 8.608357

a2
[1] 8.805906
A2
[1] 8.805906

# "simplified" expression for AIC: (AA1 - AA2) == (A1 - A2) == (a1 - a2)
AA1 <- 4 * log(s1^2) + 2 * 3
AA2 <- 4 * log(s2^2) + 2 * 4

AA1
[1] -2.743151
AA2
[1] -2.545602

• Thanks @HongOoi. I haven't understood everything concerning the regressions though. SO make a regression, meaning that we look for the function minimizing the sum of squared residuals. something like $respond = a*predictor + c$ and $respond = a*predictor + b*predictor^2 + c$ for the linear and quadratic models. Correct? Then, what is the likelhood function? – Remi.b Aug 12 '13 at 10:19