[R] Problem with MASS::fitdistr().

J C Nash pro|jcn@@h @end|ng |rom gm@||@com
Mon Apr 27 22:30:03 CEST 2020


After looking at MASS::fitdistr and fitdistrplus::fitdist, the latter seems to have
code to detect (near-)singular hessian that is almost certainly the "crash site" for
this thread. Was that package tried in this work?

I agree with Mark that writing one's own code for this is a lot of work, and I know
the folk who worked on fitdistrplus did a lot more distribution fitting problems
than I ever did, and I suspect they encountered this issue on occasions.

JN

On 2020-04-26 9:18 p.m., Mark Leeds wrote:
> it's been a looooooooong time but I vaguely remember Rvmminb computing
> gradients ( and possibly hessians )
> subject to constraints. John can say more about this but, if one is going
> to go through the anguish of
> creating a fitdstr2, then you may want to have it call Rvmminb instead of
> whatever is currently
> being called.
> 
> 
> 
> On Sun, Apr 26, 2020 at 8:55 PM Abby Spurdle <spurdle.a using gmail.com> wrote:
> 
>> I thought about this some more and realized my last suggestion is
>> unlikely to work.
>> Another possibility would be to create a new function to compute the
>> Hessian with a smaller step size, but I suspect there will be more
>> problems.
>>
>> Possibly a much simpler approach would be to:
>>
>> Modify the source for fitdistr.
>> (Copy the source and create a new function, say fitdistr2).
>>
>> Modify it not compute the Hessian in the optim call.
>> Then after the optim call, test the parameter estimates.
>> If they're very close to the boundaries (here zero), then they're
>> flagged as near-boundary cases and the fitdistr2 function returns the
>> parameter estimates without the Hessian and related info.
>> (Possibly generating a warning).
>>
>> If they're sufficiently distant, the Hessian and related info can be
>> computed in separate steps and returned.
>> (Equivalent to what it does currently).
>>
>> I note that there's at least one R package (numDeriv), and maybe more,
>> for computing the Hessian, numerically.
>>
>>
>> On Mon, Apr 27, 2020 at 9:31 AM Abby Spurdle <spurdle.a using gmail.com> wrote:
>>>
>>>> Dear Ms. Spurdle
>>>
>>> I usually refer to myself as "He".
>>> (But then, that's not the whole story...)
>>>
>>> I'm not an expert on maximum likelihood approaches.
>>> So, I apologize if the following suggestion is a poor one.
>>>
>>> Does your likelihood function have a limit, as alpha approaches zero
>> (say zero)?
>>> If so, the limit of the log-likelihood would be -Inf, would it not.
>>>
>>> You could create a function representing the likelihood or
>>> log-likelihood by wrapping your density function.
>>> The function could allow alpha/beta values equal to or below zero, and
>>> return the limit.
>>> This is mathematically incorrect, but may be sufficient for
>>> permissible estimates of the second-order partial derivatives.
>>> Depending on the shape of the likelihood function these
>>> pseudo-likelihoods maybe able to be improved...?
>>>
>>> You could then do a small modification on the source code for
>>> MASS::fitdistr, such that the user specifies the likelihood function
>>> or log-likelihood function, rather than the density...
>>>
>>> The fitdistr function is relatively complex, however, you would only
>>> need to modify a couple of lines, the lines that create the fn
>>> function...
>>
>> ______________________________________________
>> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



More information about the R-help mailing list