Coding a Gaussian log-likelihood in R

Question

I am trying to learn R by coding a Gaussian log-likelihood to solve with optim(), but after hours of sweat I am still way off the mark. (This is self-study, not homework.)

I am following the convention in many user-written functions that write a function like loglik <- function(theta, y, x) where theta is a vector of weights beta and variance sigma, y is the outcome and x is the data.

My full code with simulated data is below. Running it, you will see that my function is way off the mark compared to lm(). Can anyone give me an idea as to where I am going wrong?

# random data
set.seed(111)
y = sample(1:100,100)
x1 = sample(1:100,100)*rnorm(1,0)
x2 = sample(x1)*rnorm(1,0)
x3 = sample(x2)*rnorm(1,0)
dat = data.frame(x1,x2,x3)

# define gaussian log-likelihood
logLik <- function(theta, Y, X){
  X           <- as.matrix(X) # convert data to matrix
  k           <- ncol(X) # get the number of columns (independent vars)
  beta        <- theta[1:k] # vector of weights intialized with starting values
  expected_y  <- X %*% beta  # X is dimension (n x k) and beta is dimension (k x 1)
  sigma2      <- theta[k+1] # should be pulled from the last of the starting values vector
  LL          <- sum(dnorm(Y, mean = expected_y, sd = sigma2, log = T)) # sum of the PDF over each observation
  return(-LL)
}

Here is the output:

> optim(logLik, par=starting_values, method="Nelder-Mead", Y=y, X=dat, hessian = T)$par
[1]   0.4832514  -0.2276684  -0.3860800  32.7168490 -38.9030319
> coefficients(lm(y~x1+x2+x3))
(Intercept)          x1          x2          x3 
58.17347451 -0.06587320  0.13001865 -0.03624233

Coding a Gaussian log-likelihood in R

1 Answers