test internal math functions F32, gradF32, grad2F1, gradRegIncBeta for infinite loops

[syclik]

From @syclik on November 22, 2013 16:21

Issue #333 shows the gradregIncGamma function getting into an infinite loop. The additional functions need to be tested:

• F32
• gradF32
• grad2F1
• gradRegIncBeta

Copied from original issue: stan-dev/stan#410

[bob-carpenter]

What tests do these function need? If anyone knows, please edit the issue to add a list of tests.

[syclik]

I started working on F32, but I'm having lots of trouble. Here's the branch:
feature/issue-27-internal_math_tests

The first problem I run into is F32() isn't returning the correct results for:

  EXPECT_NEAR(0.4,
              stan::math::F32(1.0, 3.0, -1.0, 10.0, 0.5, 1.0),
              1e-6);

This is in the unit test on the branch: ./runTests.py test/unit/math/prim/scal/fun/F32_test.cpp

If there's something wrong with the input, we really need to doc the function.

[syclik]

Bump. I've gotten a first test up and it doesn't look good.

Check out branch feature/issue-27-internal_math_tests

Run:

./runTests.py test/unit/math/prim/scal/fun/F32_test.cpp

[sakrejda]

I started working on this branch.

[sakrejda]

Maybe I just haven't read the functions correctly but I was trying to find consequences to this issue in the user-facing functions and it turns out this test (that I just wrote) fails with a nan for our function.

TEST(ProbBetaBinomial, lcdf_matches_emdbook_lcdf) {
  int n = 8;
  int N = 10;
  double alpha = 3.0;
  double beta = 1.0;

  EXPECT_FLOAT_EQ(-1.849329, (stan::math::beta_binomial_lcdf(n, N, alpha, beta)));
}

With this output:

Running main() from gtest_main.cc
[==========] Running 2 tests from 1 test case.
[----------] Global test environment set-up.
[----------] 2 tests from ProbBetaBinomial
[ RUN      ] ProbBetaBinomial.cdf_log_matches_lcdf
[       OK ] ProbBetaBinomial.cdf_log_matches_lcdf (0 ms)
[ RUN      ] ProbBetaBinomial.lcdf_matches_emdbook_lcdf
test/unit/math/prim/scal/prob/beta_binomial_cdf_log_test.cpp:23: Failure
Value of: (stan::math::beta_binomial_lcdf(n, N, alpha, beta))
  Actual: nan
Expected: -1.849329
[  FAILED  ] ProbBetaBinomial.lcdf_matches_emdbook_lcdf (0 ms)
[----------] 2 tests from ProbBetaBinomial (0 ms total)

[sakrejda]

Turns out F32 also needs to terminate when the denominator for the rising factorial bits goes to 0. Once I push it will match mathematica results and have tests for that condition. Also, sign needs to be applied after exp, also fixing that. I'm not sure it was ever right except for a limited number of arguments (?)

[bob-carpenter]

Sounds good. I don't even know what this function is or where it's used, but might as well fix it if it's getting called from elsewhere. If it's not getting called anywhere, we can get rid of it unless it's there pending being exposed to the Stan language.

[sakrejda]

It gets called by the beta-binomial CDF and it's related to some functions that make my delay-time models hang so I'm just going to plow through all these if I can. I am setting up some tests against results from Mathematica so if you have any suggestions about how to do that better I appreciate suggestions. It's one of the less common instances of the generalized hypergeometric function

[sakrejda]

I think this just became a meta-issue after I dug through these and figured out what they needed. I'll open and fix smaller ones as I go.

Currently, check as fixed some issues in these functions produce wrong results or infinite loops.

• F32 assumes -1 * log(x) == log(-x), wrong result for x < 0
• F32 fix for above issue may not respect logT accumulator (but passes test?)
• F32 returns nan rather than terminating when denominator hits zero.
• grad_F32 also assumes -1 * log(x) == log(-x), wrong result for x < 0
• grad_F32 fix for above issue may not respect logT accumulator (but passes test?)
• grad_F32 can hang indefinitely
• grad_2F1 will silently return wrong results if not converging
• grad_2F1 terminates after an arbitrary number of steps
• grad_2F1 returns nan rather than terminating when denominator hits zero.
• grad_2F1 does not have a log-scale accumulator, probably should

Design across these functions:

• max_steps argument
• throw domain_error if not converging
• do ... while as in (new) F32
• return through reference to output argument (?)

Some issues are maintenance-related

• tests for values resulting in positive accumulator
• tests for values resulting in positive accumulator
• tests for values resulting in positive accumulator increment
• tests for values resulting in positive accumulator increment
• tests for throw-on-divergence, in general non-convergence should be caused by all positive a1,a2,a3,b1,b2 with large ratio or z > 1

Some issues are cosmetic but it makes it hard to figure out errors and write tests:

• grad_F32 needs doc
• grad_2F1 needs doc
• common argument naming for hypergeometric functions
• common internal naming
• common doc format
• hypergeometric functions should reference generalized hypergeometric function
• common max_steps argument, applied in same way
• prefer reference for out-type arguments over pointer (currently have both kinds!)
• grad_F32 requires a pointer-to-array[6] but doesn't doc it, only clear from reading body
• grad_reg_inc_beta needs doc
• grad_inc_beta needs doc

[hoanguc3m]

I have a problem with grad_reg_inc_beta function.
void grad_reg_inc_beta(T& g1, T& g2, const T& a, const T& b,
const T& z, const T& digammaA, const T& digammaB,
const T& digammaSum, const T& betaAB)
For example, when z = 1, I would expect that g1 is 0.
I also compare different value for z and find that it is not close to the value in the "Derivatives of the Incomplete Beta Function" (https://www.jstatsoft.org/article/view/v003i01). I dont know this function is the same with this one or not. Thanks.

[hoanguc3m]

For example, in the grad_reg_inc_beta_test.cpp

#include <stan/math/prim/scal.hpp>
#include <gtest/gtest.h>

TEST(grad_reg_inc_beta, 1) {
  double alpha = 1.0;
  double beta = 1.0;
  double y = 1.0;
  double digamma_alpha = stan::math::digamma(alpha);
  double digamma_beta = stan::math::digamma(beta);
  double digamma_sum = stan::math::digamma(alpha + beta);
  double betafunc = std::exp(stan::math::lbeta(alpha, beta));

  double g1 = 0;
  double g2 = 0;
  stan::math::grad_reg_inc_beta(g1, g2, alpha, beta, y,
                                digamma_alpha, digamma_beta, digamma_sum,
                                betafunc);
  EXPECT_FLOAT_EQ(0, g1);
  EXPECT_FLOAT_EQ(-std::numeric_limits<double>::infinity(), g2);
}

If I change alpha = 2.5, beta = 0.5. I would expect the same g1 = 0.

[sakrejda]

Thanks for pointing out the additional test case, I've been working on these but I don't think I've touched this one in any meaningful way yet. I will look at this test case but not before finishing the other issue related to these. If you're trying to use this function the best way to get help now is to ask on the users mailing list (yes, even for C++ math lib questions).

…

On Tue, Mar 28, 2017 at 5:42 AM hoanguc3m ***@***.***> wrote: For example, in the grad_reg_inc_beta_test.cpp `#include <stan/math/prim/scal.hpp> #include <gtest/gtest.h> TEST(grad_reg_inc_beta, 1) { double alpha = 1.0; double beta = 1.0; double y = 1.0; double digamma_alpha = stan::math::digamma(alpha); double digamma_beta = stan::math::digamma(beta); double digamma_sum = stan::math::digamma(alpha + beta); double betafunc = std::exp(stan::math::lbeta(alpha, beta)); double g1 = 0; double g2 = 0; stan::math::grad_reg_inc_beta(g1, g2, alpha, beta, y, digamma_alpha, digamma_beta, digamma_sum, betafunc); EXPECT_FLOAT_EQ(0, g1); EXPECT_FLOAT_EQ(-std::numeric_limits::infinity(), g2); } ` If I change alpha = 2.5, beta = 0.5. I would expect the same g1 = 0. — You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <#27 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAfA6aiPHe30jNEMT_XR4OIrYwFT8_Ozks5rqNYGgaJpZM4FTDej> .

[sakrejda]

Basically all the boxes in this issue are checked, including for grad_reg_inc_beta that relies on a hypergeometric function deep down on the inside (like most of us).