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 at 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...
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 at 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 at 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 at 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]]
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 at 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 at 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 at 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 at 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. >