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

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.

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
nu <- theta
location <- theta
scale <- theta

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  5.383121$nu
 5.779641

$location  6.259333$scale
 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  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.

• I can still estimate the variance covariance matrix by bootstrap. Yet, I'm wondering whether my estimates are valid. Apr 25, 2020 at 17:22
• I think because you're minmizing the negative log-likelihood, the standard errors are calculated by sqrt(diag(solve(fit$hessian))). As a check, I redid the maximum likelihood estimation using the bbmle package and got the same standard errors. See also here. Apr 25, 2020 at 17:45 • Thanks for spotting that, I can give you credit if you copy paste this as an answer Apr 25, 2020 at 18:48 ## 1 Answer 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
nu <- theta
location <- theta
scale <- theta

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.3530814 2.6347706 1.0937976 0.8123553


Now using bbmle:

LL2 <- function(theta){

tmp <- -sum(dpearsonIV(x, m = theta, nu = theta, location = theta, scale = theta, 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

• Is the BB optimization procedure always better than the classic optim? Apr 26, 2020 at 12:28
• Do you mean bblme and mle2? It's basically just a wrapper that makes it more convenient. Internally, it also uses optim but offers other optimizers such as nlminb, optimx etc. Apr 26, 2020 at 13:19
• Oh, I confused it with the BB optimization package. Thanks a lot! Apr 26, 2020 at 16:28