-
-
Notifications
You must be signed in to change notification settings - Fork 367
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #1433 from stan-dev/feature/issue-1426-ibeta_deriv…
…s_large_args fixes #1426. Feature/issue 1426 ibeta derivs large args
- Loading branch information
Showing
10 changed files
with
560 additions
and
100 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,94 @@ | ||
#ifndef STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDA_HPP | ||
#define STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDA_HPP | ||
|
||
#include <stan/math/prim/scal/fun/inc_beta.hpp> | ||
#include <stan/math/prim/scal/fun/inc_beta_ddb.hpp> | ||
#include <cmath> | ||
|
||
namespace stan { | ||
namespace math { | ||
|
||
template <typename T> | ||
T inc_beta_ddb(T a, T b, T z, | ||
T digamma_b, T digamma_ab); | ||
|
||
/** | ||
* Returns the partial derivative of the regularized | ||
* incomplete beta function, I_{z}(a, b) with respect to a. | ||
* The power series used to compute the deriative tends to | ||
* converge slowly when a and b are large, especially if z | ||
* approaches 1. The implementation will throw an exception | ||
* if the series have not converged within 100,000 iterations. | ||
* The current implementation has been tested for values | ||
* of a and b up to 12500 and z = 0.999. | ||
* | ||
* @tparam T scalar types of arguments | ||
* @param a a | ||
* @param b b | ||
* @param z upper bound of the integral | ||
* @param digamma_a value of digamma(a) | ||
* @param digamma_ab value of digamma(b) | ||
* @return partial derivative of the incomplete beta with respect to a | ||
* | ||
* @pre a >= 0 | ||
* @pre b >= 0 | ||
* @pre 0 <= z <= 1 | ||
*/ | ||
template <typename T> | ||
T inc_beta_dda(T a, T b, T z, | ||
T digamma_a, T digamma_ab) { | ||
using std::log; | ||
|
||
if (b > a) | ||
if ((0.1 < z && z <= 0.75 && b > 500) | ||
|| (0.01 < z && z <= 0.1 && b > 2500) | ||
|| (0.001 < z && z <= 0.01 && b > 1e5)) | ||
return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); | ||
|
||
if (z > 0.75 && a < 500) | ||
return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); | ||
if (z > 0.9 && a < 2500) | ||
return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); | ||
if (z > 0.99 && a < 1e5) | ||
return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); | ||
if (z > 0.999) | ||
return -inc_beta_ddb(b, a, 1 - z, digamma_a, digamma_ab); | ||
|
||
double threshold = 1e-10; | ||
|
||
digamma_a += 1.0 / a; // Need digamma(a + 1), not digamma(a); | ||
|
||
// Common prefactor to regularize numerator and denomentator | ||
T prefactor = (a + 1) / (a + b); | ||
prefactor = prefactor * prefactor * prefactor; | ||
|
||
T sum_numer = (digamma_ab - digamma_a) * prefactor; | ||
T sum_denom = prefactor; | ||
|
||
T summand = prefactor * z * (a + b) / (a + 1); | ||
|
||
T k = 1; | ||
digamma_ab += 1.0 / (a + b); | ||
digamma_a += 1.0 / (a + 1); | ||
|
||
while (fabs(summand) > threshold) { | ||
sum_numer += (digamma_ab - digamma_a) * summand; | ||
sum_denom += summand; | ||
|
||
summand *= (1 + (a + b) / k) * (1 + k) / (1 + (a + 1) / k); | ||
digamma_ab += 1.0 / (a + b + k); | ||
digamma_a += 1.0 / (a + 1 + k); | ||
++k; | ||
summand *= z / k; | ||
|
||
if (k > 1e5) | ||
throw std::domain_error("stan::math::inc_beta_dda did " | ||
"not converge within 100000 iterations"); | ||
} | ||
return inc_beta(a, b, z) * (log(z) + sum_numer / sum_denom); | ||
} | ||
|
||
} // math | ||
} // stan | ||
|
||
#endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,89 @@ | ||
#ifndef STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDB_HPP | ||
#define STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DDB_HPP | ||
|
||
#include <stan/math/prim/scal/fun/inc_beta.hpp> | ||
#include <stan/math/prim/scal/fun/inc_beta_dda.hpp> | ||
#include <cmath> | ||
|
||
namespace stan { | ||
namespace math { | ||
|
||
template <typename T> | ||
T inc_beta_dda(T a, T b, T z, | ||
T digamma_a, T digamma_ab); | ||
|
||
/** | ||
* Returns the partial derivative of the regularized | ||
* incomplete beta function, I_{z}(a, b) with respect to b. | ||
* The power series used to compute the deriative tends to | ||
* converge slowly when a and b are large, especailly if z | ||
* approaches 1. The implementation will throw an exception | ||
* if the series have not converged within 100,000 iterations. | ||
* The current implementation has been tested for values | ||
* of a and b up to 12500 and z = 0.999. | ||
* | ||
* @tparam T scalar types of arguments | ||
* @param a a | ||
* @param b b | ||
* @param z upper bound of the integral | ||
* @param digamma_a value of digamma(a) | ||
* @param digamma_ab value of digamma(b) | ||
* @return partial derivative of the incomplete beta with respect to b | ||
* | ||
* @pre a >= 0 | ||
* @pre b >= 0 | ||
* @pre 0 <= z <= 1 | ||
*/ | ||
template <typename T> | ||
T inc_beta_ddb(T a, T b, T z, | ||
T digamma_b, T digamma_ab) { | ||
using std::log; | ||
|
||
if (b > a) | ||
if ((0.1 < z && z <= 0.75 && b > 500) | ||
|| (0.01 < z && z <= 0.1 && b > 2500) | ||
|| (0.001 < z && z <= 0.01 && b > 1e5)) | ||
return -inc_beta_dda(b, a, 1 - z, digamma_b, digamma_ab); | ||
|
||
if ((z > 0.75 && a < 500) | ||
|| (z > 0.9 && a < 2500) | ||
|| (z > 0.99 && a < 1e5) | ||
|| (z > 0.999)) | ||
return -inc_beta_dda(b, a, 1 - z, digamma_b, digamma_ab); | ||
|
||
double threshold = 1e-10; | ||
|
||
// Common prefactor to regularize numerator and denomentator | ||
T prefactor = (a + 1) / (a + b); | ||
prefactor = prefactor * prefactor * prefactor; | ||
|
||
T sum_numer = digamma_ab * prefactor; | ||
T sum_denom = prefactor; | ||
|
||
T summand = prefactor * z * (a + b) / (a + 1); | ||
|
||
T k = 1; | ||
digamma_ab += 1.0 / (a + b); | ||
|
||
while (fabs(summand) > threshold) { | ||
sum_numer += digamma_ab * summand; | ||
sum_denom += summand; | ||
|
||
summand *= (1 + (a + b) / k) * (1 + k) / (1 + (a + 1) / k); | ||
digamma_ab += 1.0 / (a + b + k); | ||
++k; | ||
summand *= z / k; | ||
|
||
if (k > 1e5) | ||
throw std::domain_error("stan::math::inc_beta_ddb did " | ||
"not converge within 100000 iterations"); | ||
} | ||
|
||
return inc_beta(a, b, z) | ||
* (log(1 - z) - digamma_b + sum_numer / sum_denom); | ||
} | ||
|
||
} // math | ||
} // stan | ||
|
||
#endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,44 @@ | ||
#ifndef STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DERIVATIVES_HPP | ||
#define STAN_MATH_PRIM_SCAL_FUN_INC_BETA_DERIVATIVES_HPP | ||
|
||
#include <stan/math/prim/scal/fun/lgamma.hpp> | ||
#include <stan/math/prim/scal/fun/inc_beta.hpp> | ||
#include <boost/math/special_functions/beta.hpp> | ||
#include <cmath> | ||
|
||
namespace stan { | ||
namespace math { | ||
|
||
|
||
/** | ||
* Returns the partial derivative of the regularized | ||
* incomplete beta function, I_{z}(a, b) with respect to z. | ||
* | ||
* @tparam T scalar types of arguments | ||
* @param a a | ||
* @param b b | ||
* @param z upper bound of the integral | ||
* @return partial derivative of the incomplete beta with respect to z | ||
* | ||
* @pre a > 0 | ||
* @pre b > 0 | ||
* @pre 0 < z <= 1 | ||
*/ | ||
template <typename T> | ||
T inc_beta_ddz(T a, T b, T z) { | ||
using std::exp; | ||
using std::log; | ||
return exp((b - 1) * log(1 - z) + (a - 1) * log(z) | ||
+ lgamma(a + b) - lgamma(a) - lgamma(b)); | ||
} | ||
|
||
template <> | ||
double inc_beta_ddz(double a, double b, double z) { | ||
using boost::math::ibeta_derivative; | ||
return ibeta_derivative(a, b, z); | ||
} | ||
|
||
} // math | ||
} // stan | ||
|
||
#endif |
Oops, something went wrong.