Fitting three parameter log-normal distribution in R

Question

I want to fit three parameter log-normal distribution (See here for reference) in R.

My MWE is below:

set.seed(12345)
library(FAdist)
X <- rlnorm3(n=100, shape = 2, scale = 1.5, thres = 1)

# m: Location Parameter
# s: Scale Parameter
# t: Threshold Parameter
LL3 <- function(X, m, s, t)(1/((X-t)*s*(2*pi)^0.5))*exp(((-(log(X-t)-m)^2)/(2*s^2)))

library(MASS)
fitdistr(x=X, densfun=LL3, start=list(m=2, s=1.5, t=1))

But this code throws the following error message:

Error in stats::optim(x = c(30.9012208754183, 223.738029433835,
46.4287558537441,  :    non-finite finite-difference value [3] In addition: Warning message: In log(X - t) : NaNs produced

Is there any R package to fit three parameter distributions like three parameter Log-normal, Gamma, Weibull, and Log-logistic distributions?

Answer 1

In fact, it looks as though dlnorm3 (which is built into the FAdist package) already returns a zero probability when x<=thres, so plugging dlnorm3 straight into fitdistr appears to work fine:

set.seed(12345)
library(FAdist)
library(MASS)
X <- rlnorm3(n=100, shape = 2, scale = 1.5, thres = 1)
fitdistr(X,dlnorm3,start=list(shape = 2, scale = 1.5, thres = 1))

Results:

     shape        scale        thres   
  2.31116615   1.94366899   1.02798643 
 (0.18585476) (0.23426764) (0.01480906)

This does fail if we use the rllog3 function to generate values (we get much more extreme values):

Y <- rllog3(n=100, shape = 2, scale = 1.5, thres = 1)
fitdistr(Y,dlnorm3,start=list(shape = 2, scale = 1.5, thres = 1),
            method="Nelder-Mead")
## Error in stats::optim(x = c(10.1733112422871, 
##       310.508398424974, 1.08946140904075,  : 
##  non-finite finite-difference value [3]

Using debug(optim), it appears that if we switch to Nelder-Mead we can postpone the problem until the Hessian is computed.

If we use bbmle::mle2 instead we can get at least get the coefficients (with a warning that the Hessian can't be inverted ...)

library(bbmle)
mle2(Y~dlnorm3(m,s,t),
     data=data.frame(Y),
     start=list(m= 2, s = 1.5, t = 1),
         method="Nelder-Mead")


## Call:
## mle2(minuslogl = Y ~ dlnorm3(m, s, t), start = list(m = 2, s = 1.5, 
##     t = 1), method = "Nelder-Mead", data = data.frame(Y))

## Coefficients:
##        m        s        t 
## 4.227529 1.606202 1.001115 

## Log-likelihood: -440.27 
## Warning message:
## In mle2(Y ~ dlnorm3(m, s, t), data = data.frame(Y), start = list(m = 2,  :
##   couldn't invert Hessian

Answer 2

The error message indicates that there was problem in estimating the gradient of objective function. There are several reasons this could happen, but the most likely is that one of the parameters are becoming negative or causing negative values during the distribution fitting/optimization process (likely your threshold parameter becoming larger than your lognormal variable, in which case, the distribution should be 0 at those points... unfortunately, fitdistr doesn't know that).

The best way to fix this kind of problem is either to try different starting parameters or perhaps find a way to make the distribution be 0 in these cases within fitdistr.

Edit: Additionally, the code has other errors, so try, as Jealie suggested:

LL3 <- function(X, m, s, t)ifelse(t>=X,0,(1/((X-t)*s*(2*pi)^0.5))*exp(((-(log(X-t)-m)^2)/(2*s^2)))) 

and rlnorm3 rather than rllog3