Computing the gradient of the non-smooth function e.g. normal density as a function of variance parameters using R function grad

Question

I have a problem with computing a gradient of a non-smooth function using R.

To explain the core of my problem, I created a simplified problem. I have a log-likelihood trivariate normal density function where the mean is fixed to be zero and the diagonal entries of the variance are all one.

  library(mvtnorm)
  llk <- function(as,dat){

      sigma <-matrix(1,nrow=3,3)
      sigma[2,1] <- sigma[1,2] <-  as[1]
      sigma[3,1] <- sigma[1,3] <- as[2]
      sigma[3,2] <-  sigma[2,3] <- as[3]
      if (sum(eigen(sigma)$value < 0)>0) return(log(0))
      return(dmvnorm(dat,mean=c(0,0,0),sigma=sigma,log=TRUE))
      }

To keep my density to be defined for all the values of as, llk spits out -Inf if the variance is not a positive definite matrix. Now, I want to evaluate the gradient of llk with respect to the three off-diagonal entries of the variance at (0.5,0.2,0.3) and at (0.9,-0.146,0.3).

Using the grad function which is in numDeriv library, I can compute the gradient at (0.5,0.2,0.3) as:

  library(numDeriv)
  dat <- c(1,-1,2)
  as <- c(0.5,0.2,0.3)
  grad(llk,x=as,dat=dat)

 [1] -4.218858  4.533953 -6.128893

However, this method does not work for the gradient at (0.9,-0.146,0.3). To see this:

  as <- c(0.9,-0.146,0.3)
  grad(llk,x=as,dat=dat)

 > Error in grad.default(llk, x = as, dat = dat) : 
   function returns NA at 9e-051.46e-053e-05 distance from x.

I know this is because that the second variance matrix is close to non-positive definite.

Is there easy way to fix this problem? Thank you.

Computing the gradient of the non-smooth function e.g. normal density as a function of variance parameters using R function grad

0 Answers0