# 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[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.

• I can still estimate the variance covariance matrix by bootstrap. Yet, I'm wondering whether my estimates are valid. – Andrew Apr 25 '20 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. – COOLSerdash Apr 25 '20 at 17:45 • Thanks for spotting that, I can give you credit if you copy paste this as an answer – Andrew Apr 25 '20 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[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

• Is the BB optimization procedure always better than the classic optim? – Andrew Apr 26 '20 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. – COOLSerdash Apr 26 '20 at 13:19
• Oh, I confused it with the BB optimization package. Thanks a lot! – Andrew Apr 26 '20 at 16:28