2
$\begingroup$

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.

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.