Mathematica help again, @Vincent?

Developers
1

I wrote a few different calculations of d(gamma_p)/da and I’d like to document where the different versions perform best and where I should use limits/approximations. I know I have at least one where the current version does as badly as 1.0e-3 for error and the new version is much better.

I was trying to modify Vincent’s Mathematica function for computing gradients. For d(gamma_p(a,z)/da I’d like to get it to return the point the gradient was calculated at along with the gradient itself to make testing less confusing in C++. So the table I’m hoping for has the columns:

  a        z       d(gamma_p)/da   d(gamma_p)/dz
...      ...                ....             ...

Any chance I could get a hand with getting this to work? So far I’ve only managed to screw it up.

https://www.wolframcloud.com/objects/86879aa3-c1cc-4265-9ab6-407ec35446ab
2

Sorry, never mind, I got it. It’s all coming back to me now!

3

Hi Sakrejda,

I just discover your post 21h after you wrote it. Please do not hesitate to join me directly by email (picaud.vincent at gmail.com), that is safer for a quick answer. I hope you solved your problem now.

Concerning numerical computation of special functions I also would like to point this awesome blog and the associated Arb C-library which is certainly one of the best lib in this field. Surely a good source of inspiration concerning algorithms avalability and implementation.

-Vincent

4

Thanks Vincent, no worries, it was good for me to remember some Mathematica anyway. I did figure my problem out, I just needed a grid of test values for d(gamma_p(a,z))/da and I used your code along with some clunky code of mine: https://www.wolframcloud.com/objects/86879aa3-c1cc-4265-9ab6-407ec35446ab

The Arb library looks really cool for most things but I don’t see the particular derivative I’m working on in it, although I could look at how the components are implemented. I’m actually not sure how reliable the Mathematica derivatives are but they’re good enough to demonstrate the problem I’m working on: https://github.com/stan-dev/math/issues/525

Krzysztof

5

Thank you, again you are welcome

I just had a look on the Mathematica code and I noticed the way the “b” table is created before export:

b =NumberForm[Table[g[i], {i,mg}] ,18, NumberFormat → (Row[{#1, “e”, #3}] &)]

This gives something like

{0.0001e​,0.0001e​,-8.62594024578651e​}, ⋯9598⋯ ,{19.7501e​,29.7501e​,-0.01081252690344368e​}…

IMHO this can be a source of bug during parsing like “0.1e” is not a valid string to represent a float in C++.
I would suggest to replace it by

b = Map[CForm, Table[g[i], {i, mg}], {2}]

Concerning Mathematica I think it is quite reliable.

Arb is an impressive lib, implementing the latest algorithms, hence I wanted to mention it. Like it is coded in C it also allows easier integration with C++ to test some Stan code (I mean to do tests “at home”, not using Arb in the stan unit testing framework!).

To be sure, the problem is the numerical computation of d(a->GammaRegularized[a,z])?
My knowledge on the GammaRegularized function are sparse and I haven’t had the time to read the Stan code.
If I have an idea concerning that or a link to the literature I will post it here

-Vincent

6

It’s actually d(a->(1-GammaRegularized[a,z]))—derivative of the lower incomplete regularized gamma with respect to its first argument. I worked from a 1979 paper by Walter Gauschi on calculating the lower incomplete regularized gamma itself as a starting point so I’m sure we could do better (https://www.cs.purdue.edu/homes/wxg/selected_works/section_02/068.pdf). The paper has some good hints about what the major concerns are.

7

Ok. I also have found this one from Nico Temme, 1994. It also contains the Pascal source code and cites the W. Gauschi paper.

8

That’s great, the Temme paper mentions a nugget for a << 1! K

9

This raises the question: Should we be utilizing Arb in Stan Math’s unit tests?

10

I have not used Arb C lib yet.
Like it is a C lib, it would certainly allow an easy integration with Stan unit tests, however it also would add an extra dependency.
My personal option is that this does not worth it: special function computation is a very specialized and local stuff. Once it works I do not think it will be modified a lot of times and if there is a problem/bug, one can easily extract the involved code for improvements/debugging.

11

I looked at some of Arb library tests and it would be pretty easy to use (everything is really neatly wrapped up as far as I can tell. I agree with Vincent that it probably doesn’t bring too much value if we add it to math testing but as a source of values for special functions alternative to Mathematica I think it’s priceless.

For the hypergeometric functions and for the gamma derivative I recently fixed I ended up calculating values over a grid from a reference (Mathematica) to get an overview of where our function was likely failing and to choose where to branch for implementations. The one part I was really uncomfortable with was not having a good reference beyond Mathematica.

1 Like
12

Off topic, but just a few days ago I had to use Abaqus (an infinite $ finite element package) to check some calculations I was doing. It felt dirty treating some closed source package like an oracle like that haha. If it have given inconsistent results, there really wasn’t much I could do to investigate.

In the end it helped me find an inconsistency in a coordinate transformation I was doing though. Problem solved, but Pride - 1.