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.
sqrt(diag(solve(fit$hessian)))
. As a check, I redid the maximum likelihood estimation using thebbmle
package and got the same standard errors. See also here. $\endgroup$ – COOLSerdash Apr 25 '20 at 17:45