0
$\begingroup$

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:

1) The estimates between the optim and lm are different 2) 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.

EDIT

Thanks the suggestions provided by Tim and the solution from here it tunes out there were two mistakes:

1) in the llike function

 sigma_sq <- par[1+n]
 #replace by
 sigma_sq <- par[1+n]^2

2) the optimization algorithm:

Change from "BFGS" to "L-BFGS-B" and obtain the desired results.

$\endgroup$
1
$\begingroup$

Few comments:

  1. You don't need to implement the normal density by hand, you can just use build-in function. Using it is a good way to make sure that the results are not due to typo in the code.
llike <- function(par, y, X) {
    n <- ncol(X)
    beta <- par[1:n]
    sigma <- sqrt(par[1+n])
    -sum(dnorm(y, X %*% beta, sigma, log=TRUE))
}
  1. As you could see from the above code, likelihood is not distribution for the errors, but distribution for the data. In linear regression model, we predict mean of the normal distribution given some features.
  2. If you used the above code, you'd also see that the estimates are not exactly the same. Moreover, if you used different optimization algorithm, they could change again. lm function uses ordinary least squares for estimating the parameters, you used BFSG algorithm, there is no guarantee that they will give you exactly the same results.
  3. For maximizing the likelihood you need only the unnormalized log-densities, i.e. $-(x-\mu)/2\sigma^2$.

Check also the Maximum Likelihood Estimation (MLE) in layman terms thread.

$\endgroup$
3
  • $\begingroup$ thx for the reply. 1) I saw this approach, but wanted to write out the entire logl also for exercise as well. It's a important technique and you never know what densitis come your way =). 2) year, your absolutly right. 3) In fact i tried out different algorithms and got completly different results. But you dont see any obvious error in my procedure? $\endgroup$ – user2550228 Jun 3 '20 at 20:58
  • $\begingroup$ bc if no, i'll close the topic. $\endgroup$ – user2550228 Jun 3 '20 at 20:59
  • $\begingroup$ @user2550228 the general approach would be correct, but I mentioned some issues in my answer. See also the link I added to my answer for more details and friendly introduction to MLE. $\endgroup$ – Tim Jun 4 '20 at 6:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.