Edit: A link to the data and nlm() tracing information has been added, and the code has been changed in accordance with suggestions from the comments
Download the data from here. (Dropbox link).
Download the nlm() tracing information from here. (Dropbox link)
I am working on a problem to find the maximum likelihood estimators $\left( \hat\sigma_n^2,\hat l,\hat\sigma_f^2 \right)$ for some data. The log likelihood is
$$\log \mathcal{L} \left( \sigma_n^2,l,\sigma_f^2 \right)= -\frac{1}{2} y^T \left(\mathbf{K}+\sigma_n^2\mathbf{I} \right)^{-1}y – \frac{1}{2} \log \left| \mathbf{K}+\sigma_n^2\mathbf{I} \right| – \frac{n}{2} \log2\pi$$
$y$ is a $n \times 1$ vector obtained from my data (response variables). $\mathbf{I}$ is the identity matrix. $\mathbf{K}$ is a $n\times n$ matrix depending on $\sigma_f^2$ and $l$; its matrix elements are:
$$K_{pq} = \sigma_f^2 \exp\left( – \frac{1}{2l^2} \left| x_p -x_q \right|^2 \right)$$
$x_i$ are $3\times1$ vectors (independent variables).
I tried to find the maximum likelihood estimators in R by minimizing the negative log likelihood. The functions I tried using were nlm and optim. Both gave the same error – NA/Inf were produced. How do I proceed to calculate the estimators? I would prefer more accessible answers. (I am only a beginner)
I will write my code below for reference. The matrix $\mathbf{A}$ in the code (and download) is given by:
$$A_{pq} = \exp \left( – \frac{1}{2} \left| x_p -x_q \right|^2 \right) $$
load("A.Rdata") # Contains matrix A to help calculate K (in the formula above).
# A is the exponential part of the formula. (but without the 'l squared' part)
load("y.Rdata") # Contains the response variable y
num_unique <- 786
Calculate_K_plus <- function(vect){
sn2 <- vect[1]*vect[1]
exponent <- 1/(vect[2]*vect[2])
sf2 <- vect[3]*vect[3]
B <- A^exponent
B <- sf2 * B
B <- B + sn2*diag(num_unique)
return(B)}
minus_log_likelihood <- function(vect){
K_plus <- Calculate_K_plus(vect)
K_plus_inv <- solve(K_plus)
out = 0.5 * ( t(y) %*% K_plus_inv %*% y) + 0.5 * log(det(K_plus)) + (num_unique/2)*log(2*pi)
return(out)}
nlm(minus_log_likelihood,c(1,1,1))
Best Answer
This got too long for a comment. The problem with your function is that the determinant of
K_plus
was getting infinite or zero very quickly. I tweaked your function to calculate the log-determinant directly. I then usedoptim
with different methods as well asnlm
to search for the maximum likelihood estimates. The algorithms converged without problems. I also included the code to calculate the standard errors and confidence intervals based on the Hessian. All algorithms give very similar estimates.The estimates are: $\widehat{\sigma}_{n}=5.30,\hat{l}=5.12,\widehat{\sigma}_{f}=45.01$.
The code is: