Negative entries on the diagonal of the variance covariance matrix after MLE estimation of a Pearson Type 4 distribution

Question

I am trying to estimate the parameters of a Pearson Type 4 distribution using maximum likelihood. At the estimated values some of the diagonal entries of the variance-covariance matrix are not positive.

Could you please help me to solve this problem?

Take the following R example (I draw partially from the package PearsonDS):

The density is

dpearson4 <- function (x, m, nu, location, scale, log = FALSE) 
{
        k <- 2* Re(gsl::lngamma_complex(m + (nu/2) *(0+1i) ) ) - lgamma(m) - log(scale) - lgamma(m - 0.5) - lgamma(0.5)

        return(exp(k - m * log(1 + ((x - location)/scale)^2) - nu * atan((x - location)/scale)))

}

while the log-likelihood function can be written as:

LL <- function(theta, x){
  m <- theta[1]
  nu <- theta[2]
  location <- theta[3]
  scale <- theta[4]

  tmp <- -sum(log(dpearson4(x, m, nu, location, scale, log = FALSE)))

  if (is.na(tmp)) +Inf else tmp
  return(sum(tmp))
}

I generate a dummy dataset as follows (I use the rpearson function in the PearsonDS package)

set.seed(123)
x <- rpearsonIV(1000, 5, 5, 6, 6)

I start my search setting seeds using the built in ML function in PearsonDS

param <- pearsonFitML(x)[-1] 

which gives,

>param
$m
[1] 5.383121

$nu
[1] 5.779641

$location
[1] 6.259333

$scale
[1] 6.041999

However, this function does not return the Hessian so to estimate the variance-covariance matrix I run the maximum likelihood algorithm using the output of pearsonFitML as seeds

control.list <- list(maxit = 100000, factr=1e-12)#, fnscale = 1000)
  fit <- optim(par = param, 
               fn = LL, 
               hessian = TRUE, 
               method = "L-BFGS-B",
               lower = c(0.51,-Inf,-Inf,0.1),
               upper = c(Inf,Inf,Inf,Inf),
               control = control.list,
               x = x)

The output is:

solve(-fit$hessian)
                 m        nu   location      scale
m        -1.830761 -3.290973 -1.2920021 -1.0307047
nu       -3.290973 -6.941729 -2.8530122 -1.6134873
location -1.292002 -2.853012 -1.1963544 -0.6144129
scale    -1.030705 -1.613487 -0.6144129 -0.6599055
qr(fit$hessian)$rank
[1] 4
fit$par
       m       nu location    scale 
5.383123 5.779641 6.259333 6.042000 

The standard errors of the parameters are (approximately) equal to the square root of the diagonal entry of the inverse of the negative Hessian at the true value. Does this invalidate my estimate or is there an issue somewhere else? Thank you.

Answer 1

Because you're minimizing the negative log-likelihood, the standard errors are calculated by sqrt(diag(solve(fit$hessian))). Here is a short check using the bbmle package. First the original:

library(PearsonDS)
library(bbmle)

dpearson4 <- function (x, m, nu, location, scale) 
{
  k <- 2* Re(gsl::lngamma_complex(m + (nu/2) *(0+1i) ) ) - lgamma(m) - log(scale) - lgamma(m - 0.5) - lgamma(0.5)

  return(exp(k - m * log(1 + ((x - location)/scale)^2) - nu * atan((x - location)/scale)))

}

LL <- function(theta, x){
  m <- theta[1]
  nu <- theta[2]
  location <- theta[3]
  scale <- theta[4]

  tmp <- -sum(log(dpearson4(x, m, nu, location, scale)))

  if (is.na(tmp)) +Inf else tmp
  return(tmp)
}

set.seed(123)
x <- rpearsonIV(1000, 5, 5, 6, 6)

control.list <- list(maxit = 100000, factr=1e-12, fnscale = 1)

fit <- optim(par = c(2, 2, 2, 2), 
             fn = LL, 
             hessian = TRUE, 
             method = "L-BFGS-B",
             lower = c(0.51,-Inf,-Inf,0.1),
             upper = c(Inf,Inf,Inf,Inf),
             control = control.list,
             x = x)

sqrt(diag(solve(fit$hessian)))

[1] 1.3530814 2.6347706 1.0937976 0.8123553

Now using bbmle:

LL2 <- function(theta){

  tmp <- -sum(dpearsonIV(x, m = theta[1], nu = theta[2], location = theta[3], scale = theta[4], log = TRUE))
  if (is.na(tmp)) 
    +Inf
  else tmp
}

parnames(LL2) <- c("m", "nu", "location", "scale")

fit2 <- bbmle::mle2(minuslogl = LL2
            , start = c(m = 2, nu = 2, location = 2, scale = 2)
            , method = "L-BFGS-B"
            , lower = c(a = 0.5, b = -Inf, c = -Inf, d = 0.1)
            , upper = c(theta = c(Inf, Inf, Inf, Inf))
))

summary(fit2)

Coefficients:
         Estimate Std. Error z value     Pr(z)    
m         5.38541    1.35359  3.9786 6.932e-05 ***
nu        5.78422    2.63658  2.1938   0.02825 *  
location  6.26131    1.09416  5.7225 1.050e-08 ***
scale     6.04340    0.81239  7.4391 1.014e-13 ***

The standard errors are practically identical. Conveniently, the bbmle package can calculate confidence intervals based on profile likelihood:

confint(fit2)

            2.5 %    97.5 %
m        3.581845 10.276182
nu       2.680049 19.171886
location 4.695503 10.009925
scale    4.729774  8.033266