# R Different Results Calculating LLike From LM and Optim

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.

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.

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