boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Add logpdf to inverse gaussian distribution

Browse Source 
[ci skip]
This commit is contained in:
Matt Borland
2022-02-24 13:39:22 +01:00
parent 570b40e14a
commit 269bf5947f
3 changed files with 72 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
include/boost/math/distributions/inverse_gamma.hpp
View File

@@ -225,7 +225,7 @@ inline RealType logpdf(const inverse_gamma_distribution<RealType, Policy>& dist,
return 0; // random variable is infinite or so close as to make no difference.
// x * x may under or overflow, likewise our result
if (!boost::math::isfinite(x*x))
if (!(boost::math::isfinite)(x*x))
{
return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
}

37
include/boost/math/distributions/inverse_gaussian.hpp
View File

@@ -174,6 +174,43 @@ inline RealType pdf(const inverse_gaussian_distribution<RealType, Policy>& dist,
return result;
} // pdf
template <class RealType, class Policy>
inline RealType logpdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
{ // Probability Density Function
BOOST_MATH_STD_USING // for ADL of std functions
RealType scale = dist.scale();
RealType mean = dist.mean();
RealType result = 0;
static const char* function = "boost::math::logpdf(const inverse_gaussian_distribution<%1%>&, %1%)";
if(false == detail::check_scale(function, scale, &result, Policy()))
{
return result;
}
if(false == detail::check_location(function, mean, &result, Policy()))
{
return result;
}
if(false == detail::check_x_gt0(function, mean, &result, Policy()))
{
return result;
}
if(false == detail::check_positive_x(function, x, &result, Policy()))
{
return result;
}
if (x == 0)
{
return 0; // Convenient, even if not defined mathematically.
}
const RealType two_pi = boost::math::constants::two_pi<RealType>();
result = (-scale*pow(mean - x, RealType(2))/(mean*mean*x) + log(scale) - 3*log(x) - log(two_pi)) / 2;
return result;
} // pdf
template <class RealType, class Policy>
inline RealType cdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
{ // Cumulative Density Function.

35
test/test_inverse_gaussian.cpp
View File

@@ -36,6 +36,8 @@ using std::endl;
using std::setprecision;
#include <limits>
using std::numeric_limits;
#include <cmath>
using std::log;
template <class RealType>
void check_inverse_gaussian(RealType mean, RealType scale, RealType x, RealType p, RealType q, RealType tol)
@@ -207,21 +209,29 @@ BOOST_AUTO_TEST_CASE( test_main )
// formatC(SuppDists::dinverse_gaussian(1, 1, 1), digits=17) ...
BOOST_CHECK_CLOSE_FRACTION( // x = 1
pdf(w11, 1.), static_cast<double>(0.3989422804014327), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION( // x = 1
logpdf(w11, 1.), static_cast<double>(log(0.3989422804014327)), tolfeweps); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w11, 1.), static_cast<double>(0.66810200122317065), 10 * tolfeweps); // cdf
BOOST_CHECK_CLOSE_FRACTION(
pdf(w11, 0.1), static_cast<double>(0.21979480031862672), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w11, 0.1), static_cast<double>(log(0.21979480031862672)), tolfeweps); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w11, 0.1), static_cast<double>(0.0040761113207110162), 10 * tolfeweps); // cdf
BOOST_CHECK_CLOSE_FRACTION( // small x
pdf(w11, 0.01), static_cast<double>(2.0811768202028392e-19), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION( // small x
logpdf(w11, 0.01), static_cast<double>(log(2.0811768202028392e-19)), tolfeweps); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w11, 0.01), static_cast<double>(4.122313403318778e-23), 10 * tolfeweps); // cdf
BOOST_CHECK_CLOSE_FRACTION( // smaller x
pdf(w11, 0.001), static_cast<double>(2.4420044378793562e-213), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION( // smaller x
logpdf(w11, 0.001), static_cast<double>(log(2.4420044378793562e-213)), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w11, 0.001), static_cast<double>(4.8791443010851493e-219), 1000 * tolfeweps); // cdf
// 4.8791443010859224e-219 versus 4.8791443010851493e-219 so still 14 decimal digits.
@@ -240,25 +250,35 @@ BOOST_AUTO_TEST_CASE( test_main )
BOOST_CHECK_CLOSE_FRACTION(
pdf(w11, 0.5), static_cast<double>(0.87878257893544476), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w11, 0.5), static_cast<double>(log(0.87878257893544476)), tolfeweps); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w11, 0.5), static_cast<double>(0.3649755481729598), tolfeweps); // cdf
BOOST_CHECK_CLOSE_FRACTION(
pdf(w11, 2), static_cast<double>(0.10984782236693059), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w11, 2), static_cast<double>(log(0.10984782236693059)), tolfeweps); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w11, 2), static_cast<double>(.88547542598600637), tolfeweps); // cdf
BOOST_CHECK_CLOSE_FRACTION(
pdf(w11, 10), static_cast<double>(0.00021979480031862676), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w11, 10), static_cast<double>(log(0.00021979480031862676)), tolfeweps); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w11, 10), static_cast<double>(0.99964958546279115), tolfeweps); // cdf
BOOST_CHECK_CLOSE_FRACTION(
pdf(w11, 100), static_cast<double>(2.0811768202028246e-25), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w11, 100), static_cast<double>(log(2.0811768202028246e-25)), tolfeweps); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w11, 100), static_cast<double>(1), tolfeweps); // cdf
BOOST_CHECK_CLOSE_FRACTION(
pdf(w11, 1000), static_cast<double>(2.4420044378793564e-222), 10 * tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w11, 1000), static_cast<double>(log(2.4420044378793564e-222)), 10 * tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w11, 1000), static_cast<double>(1.), tolfeweps); // cdf
@@ -295,26 +315,37 @@ BOOST_AUTO_TEST_CASE( test_main )
// ===================
BOOST_CHECK_CLOSE_FRACTION( // formatC(SuppDists::dinvGauss(1, 2, 3), digits=17) "0.47490884963330904"
pdf(w23, 1.), static_cast<double>(0.47490884963330904), tolfeweps ); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w23, 1.), static_cast<double>(log(0.47490884963330904)), tolfeweps ); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
pdf(w23, 0.1), static_cast<double>(2.8854207087665401e-05), tolfeweps * 2); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w23, 0.1), static_cast<double>(log(2.8854207087665401e-05)), tolfeweps * 2); // logpdf
//2.8854207087665452e-005 2.8854207087665401e-005
BOOST_CHECK_CLOSE_FRACTION(
pdf(w23, 10.), static_cast<double>(0.0019822751498574636), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w23, 10.), static_cast<double>(log(0.0019822751498574636)), tolfeweps); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
pdf(w23, 10.), static_cast<double>(0.0019822751498574636), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w23, 10.), static_cast<double>(log(0.0019822751498574636)), tolfeweps); // logpdf
// Bigger changes in mean and scale.
inverse_gaussian w012(0.1, 2);
BOOST_CHECK_CLOSE_FRACTION(
pdf(w012, 1.), static_cast<double>(3.7460367141230404e-36), tolfeweps ); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w012, 1.), static_cast<double>(log(3.7460367141230404e-36)), tolfeweps ); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w012, 1.), static_cast<double>(1), tolfeweps ); // pdf
inverse_gaussian w0110(0.1, 10);
BOOST_CHECK_CLOSE_FRACTION(
pdf(w0110, 1.), static_cast<double>(1.6279643678071011e-176), 100 * tolfeweps ); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w0110, 1.), static_cast<double>(log(1.6279643678071011e-176)), 100 * tolfeweps ); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w0110, 1.), static_cast<double>(1), tolfeweps ); // cdf
BOOST_CHECK_CLOSE_FRACTION(
@@ -323,6 +354,8 @@ BOOST_AUTO_TEST_CASE( test_main )
BOOST_CHECK_CLOSE_FRACTION(
pdf(w0110, 0.1), static_cast<double>(39.894228040143268), tolfeweps ); // pdf
BOOST_CHECK_CLOSE_FRACTION(
logpdf(w0110, 0.1), static_cast<double>(log(39.894228040143268)), tolfeweps ); // logpdf
BOOST_CHECK_CLOSE_FRACTION(
cdf(w0110, 0.1), static_cast<double>(0.51989761564832704), 10 * tolfeweps ); // cdf