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:


# 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))

res <- optim(par=c(rep(1, 6)), llike, method="BFGS", hessian=T, y=y, X=X)

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:


Now when I calculate the logliklihood:

-llike(c(res$par[1:5], sqrt(res$par[6])/25), y, X)

compared to:

logLik(lm(y ~ X-1))

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.


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.


1 Answer 1


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.

  • $\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$ Jun 3, 2020 at 20:58
  • $\begingroup$ bc if no, i'll close the topic. $\endgroup$ Jun 3, 2020 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, 2020 at 6:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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