Bugfix/issue 487 hypergeometric functions #488

sakrejda · 2017-02-10T22:08:36Z

Submisison Checklist

Run unit tests: ./runTests.py test/unit
Run cpplint: make cpplint
Declare copyright holder and open-source license: see below

Summary:

There are currently bugs in F32 and ... I think it's grad_F32, and some potential bugs in grad_2F1 that I haven't demonstrated yet (it can terminate at an arbitrary point). See issue #487.

Intended Effect:

Fix bug in F32 and grad_F32, add some doc and tests.

How to Verify:

Run new tests versus develop to see them fail, they pass with this branch. On develop at least one of the tests hangs on my machine. Test values come from Mathematica and R's hypergeom package. The tests that fail against develop are:

cd ${STAN-MATH}
git checkout develop
git pull origin develop
git checkout bugfix/issue-487-hypergeometric-functions -- test/unit/math/prim/scal/fun/F32_test.cpp
make test/unit/math/prim/scal/fun/F32_test
./test/unit/math/prim/scal/fun/F32_test
git checkout bugfix/issue-487-hypergeometric-functions -- test/unit/math/prim/scal/prob/beta_binomial_cdf_log_test.cpp
make tests/unit/math/prim/scal/prob/beta_binomial_cdf_log_test
./test/unit/math/prim/scal/prob/beta_binomial_cdf_log_test

Side Effects:

betabinomial_cdf is now more correct and has a test that should fail with the code on develop but pass with this branch. The test seems to represent reasonable parameter values.

Documentation:

Added where it was missing.

Reviewer Suggestions:

The tests for F32 are done as far as I'm concerned and I was going to add parallel ones for grad_F32 and grad_2F1. If you want anymore let me know
The logic tracking minus signs seems bulky, suggestions appreciated but it is as hairy as it looks
F32 throws when the accumulator is infinite because that's a known non-convergence that sometimes occurs long before

Copyright and Licensing

Copyright holder: Krzysztof Sakrejda

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

…2, added tests, including some downstream for beta_binomial_cdf.

1. F32 and grad_F32 both assumed -1*log(x) == log(-x) which is not true, fixed and added a basic test. 2. From 1., the sign-flip does need to be applied after exponentiation. 3. F32 and grad_F32 both returned NaN when denominator of the increment reached zero but they should have been terminating successfully, fixed and added test comparing to Mathematica value. 4. grad_F32 could hang indefinitely, propagated Daniel's fix from F32. 5. these functions now throw when exceeding max_steps 6. grad_2F1 used to terminate after an arbitrary number of steps, now it throws if not converging after max_steps (and has new argument with default). 7. grad_2F1 would silently return wrong result of not converging, now throws.

…row on either Inf answer or max_steps.

…ator.

sakrejda · 2017-02-12T22:25:27Z

@syclik I screwed up the doxygen file in here, I'll revert it and then point you to the problem file.

sakrejda · 2017-02-13T14:32:30Z

@syclik Could you take a look at how to make this pass make test-headers? It's grad_reg_inc_beta that triggers the failure because it includes the prim version of grad_2F1 and that function is prim so it only includes the is_nan from prim (and that fails when it gets a var. If I included the rev version of is_nan in either grad_reg_inc_beta or grad_2F1 it would pass, but both of those options seem wrong since grad_2F1 calls is_nan but it's prim and grad_reg_inc_beta doesn't call is_nan ...

betanalpha · 2017-02-13T21:28:34Z

stan/math/prim/scal/fun/F32.hpp

+     *   - <code>z</code> is less than 1
+     * This function will diverge if <code>z</code> is greater than 1.
+     * When <code>z</code> is 1, which is the case for the beta binomial
+     * cdf, it is hard to determine.


Like the other hypergeometric functions, 3F2 should converge at z = 1 subject to some conditions on the other inputs. In particular, there should be a reflection symmetry that allows you to calculate near z = 1 by reflecting to 1 - z with corresponding changes to the other inputs. This is critical for ensuring robust calculations -- see, for example, https://github.com/stan-dev/math/blob/develop/stan/math/prim/scal/fun/inc_beta_dda.hpp.

👍 , I haven't fixed all the doc yet, I have a reference for when it does converge. Do you think we should be checking for convergence conditions at the front? The original function didn't do this and I was aiming for more of a bugfix here with the intention of improving the functions as a second PR. I'll either include the improvements here or make a specific second PR.

I'm leaving the reflection symmetry for a subsequent issue #511

betanalpha · 2017-02-13T21:33:41Z

stan/math/prim/scal/fun/F32.hpp


-        tNew = exp(logT);
+        if (is_nan(p) || p == 0)


I'm a little iffy about this is_nan check. When would p be NaN? Poles, such as b1 = -k, can be checked beforehand and overflow can be avoided by computing and factoring partial terms is necessary.

So if b1=-k or b2=-k this is triggered as the exit condition and I get answers within floating point precision of Mathematica as the answer. I don't think we need an alternative calculation, although maybe checking for b1=-k || b2=-k would be clearer. Suggestions?

@betanalpha many sources seem to say the function is undefined when any denominator is a negative integer, mathematica gives an answer indistinguishable from terminating (as above) at the last defined term of the sum. For this PR I'll either leave the code matching Mathematica or throw an error, I'll leave improving this situation for issue #511 as the current implementation doesn't do any better than this. Sound ok?

betanalpha · 2017-02-13T21:34:26Z

stan/math/prim/scal/fun/F32.hpp

+      bool T_is_negative = false;
+      T p = 0.0;
+      do {
+        p = (a1 + k) / (b1 + k) * (a2 + k) / (b2 + k) * (a3 + k) / (k + 1);


This line definitely needs parentheses. I much prefer the older version as it makes the possible poles much more explicit.

This is a victim of too many git commits while looking for the bug, I'll revert this line to the original which was fine.

betanalpha · 2017-02-13T21:38:21Z

stan/math/prim/scal/fun/grad_2F1.hpp

-      while (fabs(tDak * (a + (k - 1)) ) > precision || k == 0) {
-          const T r = ( (a + k) / (c + k) ) * ( (b + k) / (T)(k + 1) ) * z;
-          tDak = r * tDak * (a + (k - 1)) / (a + k);
+      do {


The do/while construction is somewhat awkward and certainly uncommon in C++. If you're worried that the while loop needs to be executed at least once then you can just do a while(true) and break on the termination condition at the end of the loop.

Seems really clear for algorithms that hit the loop at least once but I'm fine changing it.

seantalts · 2017-03-20T04:30:47Z

I'm not sure about the link to kill a build, but can you click on the build "Details" and be brought to the Jenkins page and then kill it from there?

sakrejda · 2017-03-20T14:01:07Z

@syclik I think I said this above but there's no (remaining) merge conflict. Maybe it's screwed up because I merged in develop and fixed a merge conflict with develop at some point. For the next few days I'll just leave this PR so I can address the comments already on it but then I'll do a fresh PR---it looks like when I do that as a test the new PR doesn't have any of the messy history.

syclik · 2017-03-20T18:35:45Z

Yup, the merge conflict isn't the problem now. It looks like cpplint to start.

bob-carpenter · 2017-03-20T21:24:12Z

@sakrejda This is failing due to cpplint errors.

sakrejda · 2017-03-20T22:30:02Z

Ok, it also lacks the tests I said I'd put in so it's not ready to merge anyway. Thanks.

…

On Mon, Mar 20, 2017, 5:24 PM Bob Carpenter ***@***.***> wrote: @sakrejda <https://github.com/sakrejda> This is failing due to cpplint errors. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#488 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAfA6ZGhHuCcVEQC-4O-MZP_a3Gv9Czaks5rnu58gaJpZM4L96Ge> .

…ere are multiple exit conditions.

…m:stan-dev/math into bugfix/issue-487-hypergeometric-functions Huh?

… code with trimmed down copy of 3F2 gradient code. Tests pass at 1e-8. [no ci]

…hrowing.

sakrejda · 2017-03-28T18:46:48Z

This now passes lint locally, I addressed the comments above except for the ones I punted to issue #511. The tests that call the modified functions directly or indirectly have been run AFAIK and I'm waiting for the rest to complete. The upshot was that I needed to move the gradient functions to log-scale power series as well and that lets me get EXPECT_NEAR to pass with 1e-8 error threshold while setting the precision to something like 1e-10. Previously some of the old tests were passing but with a threshold of 1e-5/6 and some of the tests I added had errors up to 1e-4. It seems like it's an open question what the default accuracy should be but now that it seems to achieve accuracy as requested I can look at some of the other gradient functions as a guide. A second round of comments would be useful at this point (feel free to wait until the tests pass) for me to figure out what needs to happen to finish this off. Keep in mind that I'd like to punt overall improvements to a separate issue and keep this to bug fixes and better defined behavior.

…diffs grad_2F1. grad_2F1 might need a non-autodiff gradient for this to be more accurate in the future.

sakrejda · 2017-03-29T13:31:10Z

@seantalts for some reason the Jenkins build fails and I can't even pull up the log, any chance you coudl look at it?

bob-carpenter · 2017-03-31T18:46:50Z

Sean's out for a couple weeks. I have no problem pulling up the log. It fails in doxygen. error report below. ``` Started by upstream project "Math Pull Request" build number 677 originally caused by: GitHub pull request #488 of commit 8aa519f, no merge conflicts. [EnvInject] - Loading node environment variables. Building remotely on gelman-group-mac (osx) in workspace /Users/Shared/Jenkins/gelman-group-mac/workspace/Math Pull Request - Doc - Doxygen /usr/bin/git rev-parse --is-inside-work-tree # timeout=10 Fetching changes from the remote Git repository /usr/bin/git config remote.origin.url https://github.com/stan-dev/math # timeout=10 Fetching upstream changes from https://github.com/stan-dev/math /usr/bin/git --version # timeout=10 /usr/bin/git fetch --tags --progress https://github.com/stan-dev/math +refs/pull/*:refs/remotes/origin/pr/* /usr/bin/git rev-parse refs/remotes/origin/pr/488/merge^{commit} # timeout=10 /usr/bin/git rev-parse refs/remotes/origin/origin/pr/488/merge^{commit} # timeout=10 Checking out Revision 2639e74 (refs/remotes/origin/pr/488/merge) /usr/bin/git config core.sparsecheckout # timeout=10 /usr/bin/git checkout -f 2639e74 /usr/bin/git rev-list d02e54e # timeout=10 [Math Pull Request - Doc - Doxygen] $ /bin/sh -xe /var/folders/8k/9dmd0c0j0h58q068bszjgs2w0000gp/T/hudson1842211594034095582.sh + git clean -d -x -f Removing doc/ Removing make/local + echo CC=clang++ + echo O=0 [Math Pull Request - Doc - Doxygen] $ /bin/sh -xe /var/folders/8k/9dmd0c0j0h58q068bszjgs2w0000gp/T/hudson903899693017782406.sh + make clean-all removing test executables rm -f -r doc/api removing dependency files rm -f rm -f lib/cvodes_2.9.0/lib/libsundials_cvodes.a lib/cvodes_2.9.0/lib/libsundials_nvecserial.a lib/cvodes_2.9.0/src/cvodes/cvodea.o lib/cvodes_2.9.0/src/cvodes/cvodea_io.o lib/cvodes_2.9.0/src/cvodes/cvodes.o lib/cvodes_2.9.0/src/cvodes/cvodes_band.o lib/cvodes_2.9.0/src/cvodes/cvodes_bandpre.o lib/cvodes_2.9.0/src/cvodes/cvodes_bbdpre.o lib/cvodes_2.9.0/src/cvodes/cvodes_dense.o lib/cvodes_2.9.0/src/cvodes/cvodes_diag.o lib/cvodes_2.9.0/src/cvodes/cvodes_direct.o lib/cvodes_2.9.0/src/cvodes/cvodes_io.o lib/cvodes_2.9.0/src/cvodes/cvodes_sparse.o lib/cvodes_2.9.0/src/cvodes/cvodes_spbcgs.o lib/cvodes_2.9.0/src/cvodes/cvodes_spgmr.o lib/cvodes_2.9.0/src/cvodes/cvodes_spils.o lib/cvodes_2.9.0/src/cvodes/cvodes_sptfqmr.o lib/cvodes_2.9.0/src/nvec_ser/nvector_serial.o lib/cvodes_2.9.0/src/sundials/sundials_band.o lib/cvodes_2.9.0/src/sundials/sundials_dense.o lib/cvodes_2.9.0/src/sundials/sundials_direct.o lib/cvodes_2.9.0/src/sundials/sundials_iterative.o lib/cvodes_2.9.0/src/sundials/sundials_math.o lib/cvodes_2.9.0/src/sundials/sundials_nvector.o lib/cvodes_2.9.0/src/sundials/sundials_pcg.o lib/cvodes_2.9.0/src/sundials/sundials_sparse.o lib/cvodes_2.9.0/src/sundials/sundials_spbcgs.o lib/cvodes_2.9.0/src/sundials/sundials_spfgmr.o lib/cvodes_2.9.0/src/sundials/sundials_spgmr.o lib/cvodes_2.9.0/src/sundials/sundials_sptfqmr.o removing generated test files rm -f + make doxygen mkdir -p doc/api doxygen doxygen/doxygen.cfg /Users/Shared/Jenkins/gelman-group-mac/workspace/Math Pull Request - Doc - Doxygen/stan/math/prim/scal/fun/grad_2F1.hpp:16: error: argument 'g' of command @param is not found in the argument list of stan::math::grad_2F1(T &g_a1, T &g_b1, const T &a1, const T &a2, const T &b1, const T &z, const T &precision=1e-10, int max_steps=1e5) (warning treated as error, aborting now) make: *** [doxygen] Error 1 Build step 'Execute shell' marked build as failure [WARNINGS] Skipping publisher since build result is FAILURE Finished: FAILURE ```

…r args.

sakrejda · 2017-04-03T23:04:05Z

I have no idea why everyone else can see the Jenkins builds on this PR and I can't I'm going to close it and start a new one.

syclik and others added 10 commits December 1, 2016 00:31

Updating F32

8ae7102

Added doc for F32; updating tests

c03e067

Added some failing tests.

c2562f1

At other branch.

bde80f2

Fixed termination condition and application of sign (after exp) on F3…

e5bb342

…2, added tests, including some downstream for beta_binomial_cdf.

Now it passes tests, still needs more tests.

209c8c7

Correctly pass the sign around, add tests for multiple sign flips, th…

7f9f1c3

…row on either Inf answer or max_steps.

Both F32 and grad_F32 use same mechanism to track sign of log-accumul…

cd3320d

…ator.

Linted.

bbe1b40

sakrejda requested review from syclik and betanalpha February 10, 2017 22:08

sakrejda self-assigned this Feb 10, 2017

sakrejda added the bug label Feb 10, 2017

sakrejda added this to the 2.14.0++ milestone Feb 10, 2017

sakrejda removed their assignment Feb 11, 2017

sakrejda added 7 commits February 11, 2017 12:44

Fix code tag.

0c9f9e9

Doc fix, checked doc compiles.

558e46a

Fixed finding is_nan.

6cac9ee

Needed to include rev version of is_nan, includes prim.

a1dd929

Fixed doc and checked doxygen.

f405369

Lint.

28a2abb

So back to this...

1d103a9

Fixed doxygen.

a809384

betanalpha reviewed Feb 13, 2017

View reviewed changes

Merge branch 'develop' into bugfix/issue-487-hypergeometric-functions

884168d

sakrejda added 2 commits March 20, 2017 10:02

Whut?

9935b09

Added tests and fixed convergence check functions for 3F2 and 2F1.

7c9fa78

sakrejda mentioned this pull request Mar 20, 2017

Hypergeometric function improvements #511

Closed

5 tasks

Old calculation makes poles clearer.

b12a563

sakrejda added 7 commits March 28, 2017 11:48

Moved grad_F32 over to log-based calculations.

034b558

Reverted precision request to default and used while(true) loop as th…

59dced8

…ere are multiple exit conditions.

Merge branch 'bugfix/issue-487-hypergeometric-functions' of github.co…

3e8987e

…m:stan-dev/math into bugfix/issue-487-hypergeometric-functions Huh?

Up requested accuracy in tests so current code fails, replace current…

e3d1fe6

… code with trimmed down copy of 3F2 gradient code. Tests pass at 1e-8. [no ci]

Exercise more of the branches in grad_2F1.hpp for sign-flipping and t…

420934c

…hrowing.

Whitespace lint.

c8694f1

Lint.

c698fa3

Adjusted tolerance 1.05e-6 vs. 1e-6 on 3rd order diff test that auto-…

8aa519f

…diffs grad_2F1. grad_2F1 might need a non-autodiff gradient for this to be more accurate in the future.

sakrejda added 3 commits April 3, 2017 16:55

is_polynomial at zero numerator args, is_undefined at zero denominato…

657b936

…r args.

Pass cpplint.

0722f30

Dropped extra iostream includes and passed all tests in test/unit/math.

82aec7a

sakrejda closed this Apr 3, 2017

seantalts modified the milestones: 2.14.0++, v2.15.0 Apr 14, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Bugfix/issue 487 hypergeometric functions #488

Bugfix/issue 487 hypergeometric functions #488

sakrejda commented Feb 10, 2017 •

edited

Loading

sakrejda commented Feb 12, 2017

sakrejda commented Feb 13, 2017

betanalpha Feb 13, 2017

sakrejda Feb 14, 2017

sakrejda Mar 20, 2017

betanalpha Feb 13, 2017

sakrejda Feb 14, 2017

sakrejda Mar 20, 2017

betanalpha Feb 13, 2017

sakrejda Feb 14, 2017

betanalpha Feb 13, 2017

sakrejda Feb 14, 2017

seantalts commented Mar 20, 2017

sakrejda commented Mar 20, 2017

syclik commented Mar 20, 2017

bob-carpenter commented Mar 20, 2017

sakrejda commented Mar 20, 2017 via email

sakrejda commented Mar 28, 2017

sakrejda commented Mar 29, 2017

bob-carpenter commented Mar 31, 2017 via email

sakrejda commented Apr 3, 2017

Bugfix/issue 487 hypergeometric functions #488

Bugfix/issue 487 hypergeometric functions #488

Conversation

sakrejda commented Feb 10, 2017 • edited Loading

Submisison Checklist

Summary:

Intended Effect:

How to Verify:

Side Effects:

Documentation:

Reviewer Suggestions:

Copyright and Licensing

sakrejda commented Feb 12, 2017

sakrejda commented Feb 13, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seantalts commented Mar 20, 2017

sakrejda commented Mar 20, 2017

syclik commented Mar 20, 2017

bob-carpenter commented Mar 20, 2017

sakrejda commented Mar 20, 2017 via email

sakrejda commented Mar 28, 2017

sakrejda commented Mar 29, 2017

bob-carpenter commented Mar 31, 2017 via email

sakrejda commented Apr 3, 2017

sakrejda commented Feb 10, 2017 •

edited

Loading