I am trying to implement a multivariate linear regression model by ML estimation, however I ran into some discrepancies.
An assumption of the model is that the residual is a standard normal variable.
I followed this(page 6) approach, but the procedure should be clear anyway:
$$ \text{pdf}_X (x,\mu,\sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} $$
$$ \text{pdf}_X (x_1,...,x_n,\mu,\sigma) = \prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x_i-\mu)^2}{2\sigma^2}} = (\frac{1}{\sqrt{2\pi\sigma^2}})^ne^{-\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}} $$
Taking the log gives us the log likelihood:
$$ log((\frac{1}{\sqrt{2\pi\sigma^2}})^ne^{-\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}}) = -\frac{n}{2}log(2\pi) -\frac{n}{2}log(\sigma^2) - \frac{\sum_{i=1}^n(y - X\beta)^2}{2\sigma^2} $$
Translating this into R Code:
set.seed(1)
# generate random data
X <- cbind(rep(1, 100/4), matrix(rnorm(100), ncol=4, byrow=T))
#true beta
beta <- c(25, 1,3,5,7)
# linear model
y <- X%*%beta + rnorm(25)
llike <- function(par, y, X){
m <- nrow(X)
n <- ncol(X)
beta <- par[1:n]
sigma_sq <- par[1+n]
e <- y - X%*%beta
loglike <- -(m/2)*log(2*pi) - (m/2)*log(sigma_sq) - ((t(e)%*%e) / (2*sigma_sq))
return(-loglike)
}
res <- optim(par=c(rep(1, 6)), llike, method="BFGS", hessian=T, y=y, X=X)
res$par
24.189 1.209 3.729 5.489 7.102 8788.176
The last term is the variance it seems excessive, but keep in mind that the "rnorm" function generates normal distributed data, with mean = 0 and standard deviation = 1. Furthermore we have 25 of those so the estimated standard deviation would actually be:
sqrt(res$par[6])/25
3.749
Now when I calculate the logliklihood:
-llike(c(res$par[1:5], sqrt(res$par[6])/25), y, X)
-45.56091
compared to:
logLik(lm(y ~ X-1))
-42.238
So there are two issues that I have:
- The estimates between the optim and lm are different
- the loglikelihood is differes.
Am I missing something with the implementation? I really need to know if there is a mistake, since this model is the bases for Lasso and Ridge Regression which are supposed to follow.