Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add trigamma and invdigamma functions #3226

Closed
wants to merge 3 commits into from
Closed

Add trigamma and invdigamma functions #3226

wants to merge 3 commits into from

Conversation

johnmyleswhite
Copy link
Member

This pull request adds the trigamma and invdigamma functions, which are the major missing instances of the psigamma function.

The tests verify that the functions return values that are roughly comparable with the values provided by R.

@johnmyleswhite
Copy link
Member Author

It would be good to get @andreasnoackjensen's feedback on this.

@JeffBezanson
Copy link
Sponsor Member

Awesome, thanks!

@andreasnoack
Copy link
Member

@johnmyleswhite I can confirm the slight difference between as121 and R/MATLAB/Mathematica. R cites: Amos, D. E. (1983). A portable Fortran subroutine for derivatives of the psi function, Algorithm 610, and it is part of SLATEC, http://www.netlib.org/slatec/src/dpsifn.f and therefore in the public domain. I tried to translate that to Julia at it seem to work. I think it is the more precise of the two algorithms and at the same time it gives all the derivatives so I suggest we use Amos' algorithm.

I don't have much to say about the inverse digamma function. I haven't encountered the function before and it also seems a little rare with around 240 hits on Google, but the implementation seems to work fine.

@johnmyleswhite
Copy link
Member Author

Do you have a working implementation of Amos' algorithm? If so, let's definitely use that.

The inverse digamma function comes up in fixed point algorithms involving the digamma function, e.g. MLE for the Dirichlet distribution.

@andreasnoack
Copy link
Member

I think it works. I'll make a pull request in a minute. There are still some open questions regarding the code, but it is easier to discuss when you can see the code.

@stevengj
Copy link
Member

These algorithms (especially invdigamma) look insanely slow compared to just chopping up the domain into a few pieces and using a minimax rational approximant on each piece. (Nowadays, Maple or Mathematica can construct minimax approximants for you, to arbitrary precision, so implementing this strategy is not so hard.)

But I supposed they can always be replaced by an optimized implementation later.

@johnmyleswhite
Copy link
Member Author

Closing this. I'll add invdigamma on its own now that we have polygamma.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants