From 8752b55fabacc9bdf2ece5dfb89a1426535dcb62 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Sun, 27 Oct 2013 03:04:36 -0400 Subject: [PATCH 1/4] added two tests that fail in gamma_q chain calculations --- src/test/agrad/rev/gamma_q_test.cpp | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/test/agrad/rev/gamma_q_test.cpp b/src/test/agrad/rev/gamma_q_test.cpp index a0a63dda69c..09e43daea7e 100644 --- a/src/test/agrad/rev/gamma_q_test.cpp +++ b/src/test/agrad/rev/gamma_q_test.cpp @@ -21,6 +21,26 @@ TEST(AgradRev,gamma_q_var_var) { b = -1.0; EXPECT_THROW(gamma_q(a,b), std::domain_error); } +TEST(AgradRevGammaQ,infLoopInVersion2_0_1_var_var) { + // FIXME: causes infinite loop in 2.0.1 gradient calcs + AVAR a = 8.01006; + AVAR b = 2.47579e+215; + AVEC x = createAVEC(a,b); + + AVAR f = gamma_q(a,b); + VEC g; + EXPECT_THROW(f.grad(x,g), std::domain_error); +} +TEST(AgradRevGammaQ,infLoopInVersion2_0_1_var_double) { + // FIXME: causes infinite loop in 2.0.1 gradient calcs + AVAR a = 8.01006; + double b = 2.47579e+215; + AVEC x = createAVEC(a); + + AVAR f = gamma_q(a,b); + VEC g; + EXPECT_THROW(f.grad(x,g), std::domain_error); +} TEST(AgradRev,gamma_q_double_var) { double a = 0.5; AVAR b = 1.0; From 5f521435524f11c1e9ff3cf892296ce83af81a79 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Thu, 21 Nov 2013 23:03:19 -0500 Subject: [PATCH 2/4] refactored code. The same implementation was in src/stan/prob/internal_math.hpp. This doesn't solve the infinite loop problem. --- src/stan/agrad/rev/gamma_q.hpp | 58 ++++++++-------------------------- 1 file changed, 14 insertions(+), 44 deletions(-) diff --git a/src/stan/agrad/rev/gamma_q.hpp b/src/stan/agrad/rev/gamma_q.hpp index b1e69418781..ba0ca04cdc0 100644 --- a/src/stan/agrad/rev/gamma_q.hpp +++ b/src/stan/agrad/rev/gamma_q.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -21,31 +22,15 @@ namespace stan { avi,bvi) { } void chain() { - - double u = stan::math::gamma_q(avi_->val_, bvi_->val_); - - double S = 0; - double s = 1; - double l = std::log(bvi_->val_); - double g = boost::math::tgamma(avi_->val_); - double dig = boost::math::digamma(avi_->val_); - - int k = 0; - double delta = s / (avi_->val_ * avi_->val_); - - while (std::fabs(delta) > 1e-6) { - S += delta; - ++k; - s *= - bvi_->val_ / k; - delta = s / ((k + avi_->val_) * (k + avi_->val_)); - } - - - avi_->adj_ += adj_ * ((1.0 - u) * ( dig - l ) + std::exp( avi_->val_ * l ) * S / g); - bvi_->adj_ -= adj_ * (std::exp(-bvi_->val_) * std::pow(bvi_->val_, avi_->val_ - 1.0) / g); + avi_->adj_ += adj_ + * stan::math::gradRegIncGamma(avi_->val_, bvi_->val_, + boost::math::tgamma(avi_->val_), + boost::math::digamma(avi_->val_)); + bvi_->adj_ -= adj_ + * boost::math::gamma_p_derivative(avi_->val_, bvi_->val_); } }; - + class gamma_q_vd_vari : public op_vd_vari { public: gamma_q_vd_vari(vari* avi, double b) : @@ -53,26 +38,10 @@ namespace stan { avi,b) { } void chain() { - - double u = stan::math::gamma_q(avi_->val_, bd_); - - double S = 0; - double s = 1; - double l = std::log(bd_); - double g = boost::math::tgamma(avi_->val_); - double dig = boost::math::digamma(avi_->val_); - - int k = 0; - double delta = s / (avi_->val_ * avi_->val_); - - while (std::fabs(delta) > 1e-6) { - S += delta; - ++k; - s *= - bd_ / k; - delta = s / ((k + avi_->val_) * (k + avi_->val_)); - } - - avi_->adj_ += adj_ * ((1.0 - u) * ( dig - l ) + std::exp( avi_->val_ * l ) * S / g); + avi_->adj_ += adj_ + * stan::math::gradRegIncGamma(avi_->val_, bd_, + boost::math::tgamma(avi_->val_), + boost::math::digamma(avi_->val_)); } }; @@ -83,7 +52,8 @@ namespace stan { a,bvi) { } void chain() { - bvi_->adj_ -= adj_ * (std::exp(-bvi_->val_) * std::pow(bvi_->val_, ad_ - 1.0) / boost::math::tgamma(ad_)); + bvi_->adj_ -= adj_ + * boost::math::gamma_p_derivative(ad_, bvi_->val_); } }; } From 43d8a8f5cc72d646e32886706512ea2e06fa702a Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Tue, 26 Nov 2013 20:56:23 -0500 Subject: [PATCH 3/4] updating stan/prob/internal_math function to throw an exception when the sum diverges --- src/stan/prob/internal_math.hpp | 54 ++++++++++++++-------------- src/test/prob/internal_math_test.cpp | 21 +++++++++++ 2 files changed, 49 insertions(+), 26 deletions(-) create mode 100644 src/test/prob/internal_math_test.cpp diff --git a/src/stan/prob/internal_math.hpp b/src/stan/prob/internal_math.hpp index 5315eaa6a8f..9cb298768d4 100644 --- a/src/stan/prob/internal_math.hpp +++ b/src/stan/prob/internal_math.hpp @@ -5,6 +5,8 @@ #include #include +#include + namespace stan { namespace math { @@ -173,33 +175,33 @@ namespace stan { } - // Gradient of the regularized incomplete gamma functions igamma(a, g) - double gradRegIncGamma(double a, double z, double g, double dig, double precision = 1e-6) - { - - using boost::math::gamma_p; - - double S = 0; - double s = 1; - double l = std::log(z); - - int k = 0; - double delta = s / (a * a); - - while (fabs(delta) > precision) - { - S += delta; - ++k; - s *= - z / k; - delta = s / ((k + a) * (k + a)); - } - - // Precomputed values - // dig -> digamma(a) - // g -> g(a) - return gamma_p(a, z) * ( dig - l ) + std::exp( a * l ) * S / g; - + // Gradient of the regularized incomplete gamma functions igamma(a, g) + // Precomputed values + // g = boost::math::tgamma(a) + // dig = boost::math::digamma(a) + double gradRegIncGamma(double a, double z, double g, double dig, + double precision = 1e-6) { + using boost::math::gamma_p; + + double S = 0; + double s = 1; + double l = std::log(z); + + int k = 0; + double delta = s / (a * a); + double last_delta; + + while (fabs(delta) > precision) { + S += delta; + last_delta = delta; + ++k; + s *= - z / k; + delta = s / ((k + a) * (k + a)); + if (fabs(delta) > fabs(last_delta)) + throw std::domain_error("stan::math::gradRegIncGamma not converging"); } + return gamma_p(a, z) * ( dig - l ) + std::exp( a * l ) * S / g; + } } diff --git a/src/test/prob/internal_math_test.cpp b/src/test/prob/internal_math_test.cpp new file mode 100644 index 00000000000..2067b584759 --- /dev/null +++ b/src/test/prob/internal_math_test.cpp @@ -0,0 +1,21 @@ +#include +#include + +TEST(ProbInternalMath, gradRegIncGamma_typical) { + double a = 0.5; + double b = 1.0; + double g = 1.77245; + double dig = -1.96351; + + EXPECT_FLOAT_EQ(0.38984156, stan::math::gradRegIncGamma(a, b, g, dig)); +} + +TEST(ProbInternalMath, gradRegIncGamma_infLoopInVersion2_0_1) { + double a = 8.01006; + double b = 2.47579e+215; + double g = 5143.28; + double dig = 2.01698; + + EXPECT_THROW(stan::math::gradRegIncGamma(a, b, g, dig), + std::domain_error); +} From 368c02ef676a6e8b9a6560fd6d1e79be3bc7d622 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Wed, 27 Nov 2013 15:01:08 -0500 Subject: [PATCH 4/4] changing error condition to when there is actually an error --- src/stan/prob/internal_math.hpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/stan/prob/internal_math.hpp b/src/stan/prob/internal_math.hpp index 9cb298768d4..7ccb5e223b5 100644 --- a/src/stan/prob/internal_math.hpp +++ b/src/stan/prob/internal_math.hpp @@ -4,8 +4,7 @@ #include #include #include - -#include +#include namespace stan { @@ -189,15 +188,13 @@ namespace stan { int k = 0; double delta = s / (a * a); - double last_delta; while (fabs(delta) > precision) { S += delta; - last_delta = delta; ++k; s *= - z / k; delta = s / ((k + a) * (k + a)); - if (fabs(delta) > fabs(last_delta)) + if (boost::math::isinf(delta)) throw std::domain_error("stan::math::gradRegIncGamma not converging"); } return gamma_p(a, z) * ( dig - l ) + std::exp( a * l ) * S / g;