boost/math/distributions/inverse_gaussian.hpp

   1 //  Copyright John Maddock 2010.
   2 //  Copyright Paul A. Bristow 2010.
   3
   4 //  Use, modification and distribution are subject to the
   5 //  Boost Software License, Version 1.0. (See accompanying file
   6 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #ifndef BOOST_STATS_INVERSE_GAUSSIAN_HPP
   9 #define BOOST_STATS_INVERSE_GAUSSIAN_HPP
  10
  11 #ifdef _MSC_VER
  12 #pragma warning(disable: 4512) // assignment operator could not be generated
  13 #endif
  14
  15 // http://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution
  16 // http://mathworld.wolfram.com/InverseGaussianDistribution.html
  17
  18 // The normal-inverse Gaussian distribution
  19 // also called the Wald distribution (some sources limit this to when mean = 1).
  20
  21 // It is the continuous probability distribution
  22 // that is defined as the normal variance-mean mixture where the mixing density is the
  23 // inverse Gaussian distribution. The tails of the distribution decrease more slowly
  24 // than the normal distribution. It is therefore suitable to model phenomena
  25 // where numerically large values are more probable than is the case for the normal distribution.
  26
  27 // The Inverse Gaussian distribution was first studied in relationship to Brownian motion.
  28 // In 1956 M.C.K. Tweedie used the name 'Inverse Gaussian' because there is an inverse
  29 // relationship between the time to cover a unit distance and distance covered in unit time.
  30
  31 // Examples are returns from financial assets and turbulent wind speeds.
  32 // The normal-inverse Gaussian distributions form
  33 // a subclass of the generalised hyperbolic distributions.
  34
  35 // See also
  36
  37 // http://en.wikipedia.org/wiki/Normal_distribution
  38 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
  39 // Also:
  40 // Weisstein, Eric W. "Normal Distribution."
  41 // From MathWorld--A Wolfram Web Resource.
  42 // http://mathworld.wolfram.com/NormalDistribution.html
  43
  44 // http://www.jstatsoft.org/v26/i04/paper General class of inverse Gaussian distributions.
  45 // ig package - withdrawn but at http://cran.r-project.org/src/contrib/Archive/ig/
  46
  47 // http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/SuppDists/html/inverse_gaussian.html
  48 // R package for dinverse_gaussian, ...
  49
  50 // http://www.statsci.org/s/inverse_gaussian.s  and http://www.statsci.org/s/inverse_gaussian.html
  51
  52 //#include <boost/math/distributions/fwd.hpp>
  53 #include <boost/math/special_functions/erf.hpp> // for erf/erfc.
  54 #include <boost/math/distributions/complement.hpp>
  55 #include <boost/math/distributions/detail/common_error_handling.hpp>
  56 #include <boost/math/distributions/normal.hpp>
  57 #include <boost/math/distributions/gamma.hpp> // for gamma function
  58 // using boost::math::gamma_p;
  59
  60 #include <boost/math/tools/tuple.hpp>
  61 //using std::tr1::tuple;
  62 //using std::tr1::make_tuple;
  63 #include <boost/math/tools/roots.hpp>
  64 //using boost::math::tools::newton_raphson_iterate;
  65
  66 #include <utility>
  67
  68 namespace boost{ namespace math{
  69
  70 template <class RealType = double, class Policy = policies::policy<> >
  71 class inverse_gaussian_distribution
  72 {
  73 public:
  74    typedef RealType value_type;
  75    typedef Policy policy_type;
  76
  77    inverse_gaussian_distribution(RealType mean = 1, RealType scale = 1)
  78       : m_mean(mean), m_scale(scale)
  79    { // Default is a 1,1 inverse_gaussian distribution.
  80      static const char* function = "boost::math::inverse_gaussian_distribution<%1%>::inverse_gaussian_distribution";
  81
  82      RealType result;
  83      detail::check_scale(function, scale, &result, Policy());
  84      detail::check_location(function, mean, &result, Policy());
  85    }
  86
  87    RealType mean()const
  88    { // alias for location.
  89       return m_mean; // aka mu
  90    }
  91
  92    // Synonyms, provided to allow generic use of find_location and find_scale.
  93    RealType location()const
  94    { // location, aka mu.
  95       return m_mean;
  96    }
  97    RealType scale()const
  98    { // scale, aka lambda.
  99       return m_scale;
 100    }
 101
 102    RealType shape()const
 103    { // shape, aka phi = lambda/mu.
 104       return m_scale / m_mean;
 105    }
 106
 107 private:
 108    //
 109    // Data members:
 110    //
 111    RealType m_mean;  // distribution mean or location, aka mu.
 112    RealType m_scale;    // distribution standard deviation or scale, aka lambda.
 113 }; // class normal_distribution
 114
 115 typedef inverse_gaussian_distribution<double> inverse_gaussian;
 116
 117 template <class RealType, class Policy>
 118 inline const std::pair<RealType, RealType> range(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)
 119 { // Range of permissible values for random variable x, zero to max.
 120    using boost::math::tools::max_value;
 121    return std::pair<RealType, RealType>(static_cast<RealType>(0.), max_value<RealType>()); // - to + max value.
 122 }
 123
 124 template <class RealType, class Policy>
 125 inline const std::pair<RealType, RealType> support(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)
 126 { // Range of supported values for random variable x, zero to max.
 127   // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 128    using boost::math::tools::max_value;
 129    return std::pair<RealType, RealType>(static_cast<RealType>(0.),  max_value<RealType>()); // - to + max value.
 130 }
 131
 132 template <class RealType, class Policy>
 133 inline RealType pdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
 134 { // Probability Density Function
 135    BOOST_MATH_STD_USING  // for ADL of std functions
 136
 137    RealType scale = dist.scale();
 138    RealType mean = dist.mean();
 139    RealType result = 0;
 140    static const char* function = "boost::math::pdf(const inverse_gaussian_distribution<%1%>&, %1%)";
 141    if(false == detail::check_scale(function, scale, &result, Policy()))
 142    {
 143       return result;
 144    }
 145    if(false == detail::check_location(function, mean, &result, Policy()))
 146    {
 147       return result;
 148    }
 149    if(false == detail::check_positive_x(function, x, &result, Policy()))
 150    {
 151       return result;
 152    }
 153
 154    if (x == 0)
 155    {
 156      return 0; // Convenient, even if not defined mathematically.
 157    }
 158
 159    result =
 160      sqrt(scale / (constants::two_pi<RealType>() * x * x * x))
 161     * exp(-scale * (x - mean) * (x - mean) / (2 * x * mean * mean));
 162    return result;
 163 } // pdf
 164
 165 template <class RealType, class Policy>
 166 inline RealType cdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
 167 { // Cumulative Density Function.
 168    BOOST_MATH_STD_USING  // for ADL of std functions.
 169
 170    RealType scale = dist.scale();
 171    RealType mean = dist.mean();
 172    static const char* function = "boost::math::cdf(const inverse_gaussian_distribution<%1%>&, %1%)";
 173    RealType result = 0;
 174    if(false == detail::check_scale(function, scale, &result, Policy()))
 175    {
 176       return result;
 177    }
 178    if(false == detail::check_location(function, mean, &result, Policy()))
 179    {
 180       return result;
 181    }
 182    if(false == detail::check_positive_x(function, x, &result, Policy()))
 183    {
 184      return result;
 185    }
 186    if (x == 0)
 187    {
 188      return 0; // Convenient, even if not defined mathematically.
 189    }
 190    // Problem with this formula for large scale > 1000 or small x,
 191    //result = 0.5 * (erf(sqrt(scale / x) * ((x / mean) - 1) / constants::root_two<RealType>(), Policy()) + 1)
 192    //  + exp(2 * scale / mean) / 2
 193    //  * (1 - erf(sqrt(scale / x) * (x / mean + 1) / constants::root_two<RealType>(), Policy()));
 194    // so use normal distribution version:
 195    // Wikipedia CDF equation http://en.wikipedia.org/wiki/Inverse_Gaussian_distribution.
 196
 197    normal_distribution<RealType> n01;
 198
 199    RealType n0 = sqrt(scale / x);
 200    n0 *= ((x / mean) -1);
 201    RealType n1 = cdf(n01, n0);
 202    RealType expfactor = exp(2 * scale / mean);
 203    RealType n3 = - sqrt(scale / x);
 204    n3 *= (x / mean) + 1;
 205    RealType n4 = cdf(n01, n3);
 206    result = n1 + expfactor * n4;
 207    return result;
 208 } // cdf
 209
 210 template <class RealType>
 211 struct inverse_gaussian_quantile_functor
 212 {
 213
 214   inverse_gaussian_quantile_functor(const boost::math::inverse_gaussian_distribution<RealType> dist, RealType const& p)
 215     : distribution(dist), prob(p)
 216   {
 217   }
 218   boost::math::tuple<RealType, RealType> operator()(RealType const& x)
 219   {
 220     RealType c = cdf(distribution, x);
 221     RealType fx = c - prob;  // Difference cdf - value - to minimize.
 222     RealType dx = pdf(distribution, x); // pdf is 1st derivative.
 223     // return both function evaluation difference f(x) and 1st derivative f'(x).
 224     return boost::math::make_tuple(fx, dx);
 225   }
 226   private:
 227   const boost::math::inverse_gaussian_distribution<RealType> distribution;
 228   RealType prob;
 229 };
 230
 231 template <class RealType>
 232 struct inverse_gaussian_quantile_complement_functor
 233 {
 234     inverse_gaussian_quantile_complement_functor(const boost::math::inverse_gaussian_distribution<RealType> dist, RealType const& p)
 235     : distribution(dist), prob(p)
 236   {
 237   }
 238   boost::math::tuple<RealType, RealType> operator()(RealType const& x)
 239   {
 240     RealType c = cdf(complement(distribution, x));
 241     RealType fx = c - prob;  // Difference cdf - value - to minimize.
 242     RealType dx = -pdf(distribution, x); // pdf is 1st derivative.
 243     // return both function evaluation difference f(x) and 1st derivative f'(x).
 244     //return std::tr1::make_tuple(fx, dx); if available.
 245     return boost::math::make_tuple(fx, dx);
 246   }
 247   private:
 248   const boost::math::inverse_gaussian_distribution<RealType> distribution;
 249   RealType prob;
 250 };
 251
 252 namespace detail
 253 {
 254   template <class RealType>
 255   inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1)
 256   { // guess at random variate value x for inverse gaussian quantile.
 257       BOOST_MATH_STD_USING
 258       using boost::math::policies::policy;
 259       // Error type.
 260       using boost::math::policies::overflow_error;
 261       // Action.
 262       using boost::math::policies::ignore_error;
 263
 264       typedef policy<
 265         overflow_error<ignore_error> // Ignore overflow (return infinity)
 266       > no_overthrow_policy;
 267
 268     RealType x; // result is guess at random variate value x.
 269     RealType phi = lambda / mu;
 270     if (phi > 2.)
 271     { // Big phi, so starting to look like normal Gaussian distribution.
 272       //    x=(qnorm(p,0,1,true,false) - 0.5 * sqrt(mu/lambda)) / sqrt(lambda/mu);
 273       // Whitmore, G.A. and Yalovsky, M.
 274       // A normalising logarithmic transformation for inverse Gaussian random variables,
 275       // Technometrics 20-2, 207-208 (1978), but using expression from
 276       // V Seshadri, Inverse Gaussian distribution (1998) ISBN 0387 98618 9, page 6.
 277
 278       normal_distribution<RealType, no_overthrow_policy> n01;
 279       x = mu * exp(quantile(n01, p) / sqrt(phi) - 1/(2 * phi));
 280      }
 281     else
 282     { // phi < 2 so much less symmetrical with long tail,
 283       // so use gamma distribution as an approximation.
 284       using boost::math::gamma_distribution;
 285
 286       // Define the distribution, using gamma_nooverflow:
 287       typedef gamma_distribution<RealType, no_overthrow_policy> gamma_nooverflow;
 288
 289       gamma_distribution<RealType, no_overthrow_policy> g(static_cast<RealType>(0.5), static_cast<RealType>(1.));
 290
 291       // gamma_nooverflow g(static_cast<RealType>(0.5), static_cast<RealType>(1.));
 292       // R qgamma(0.2, 0.5, 1)  0.0320923
 293       RealType qg = quantile(complement(g, p));
 294       //RealType qg1 = qgamma(1.- p, 0.5, 1.0, true, false);
 295       x = lambda / (qg * 2);
 296       //
 297       if (x > mu/2) // x > mu /2?
 298       { // x too large for the gamma approximation to work well.
 299         //x = qgamma(p, 0.5, 1.0); // qgamma(0.270614, 0.5, 1) = 0.05983807
 300         RealType q = quantile(g, p);
 301        // x = mu * exp(q * static_cast<RealType>(0.1));  // Said to improve at high p
 302        // x = mu * x;  // Improves at high p?
 303         x = mu * exp(q / sqrt(phi) - 1/(2 * phi));
 304       }
 305     }
 306     return x;
 307   }  // guess_ig
 308 } // namespace detail
 309
 310 template <class RealType, class Policy>
 311 inline RealType quantile(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& p)
 312 {
 313    BOOST_MATH_STD_USING  // for ADL of std functions.
 314    // No closed form exists so guess and use Newton Raphson iteration.
 315
 316    RealType mean = dist.mean();
 317    RealType scale = dist.scale();
 318    static const char* function = "boost::math::quantile(const inverse_gaussian_distribution<%1%>&, %1%)";
 319
 320    RealType result = 0;
 321    if(false == detail::check_scale(function, scale, &result, Policy()))
 322       return result;
 323    if(false == detail::check_location(function, mean, &result, Policy()))
 324       return result;
 325    if(false == detail::check_probability(function, p, &result, Policy()))
 326       return result;
 327    if (p == 0)
 328    {
 329      return 0; // Convenient, even if not defined mathematically?
 330    }
 331    if (p == 1)
 332    { // overflow
 333       result = policies::raise_overflow_error<RealType>(function,
 334         "probability parameter is 1, but must be < 1!", Policy());
 335       return result; // std::numeric_limits<RealType>::infinity();
 336    }
 337
 338   RealType guess = detail::guess_ig(p, dist.mean(), dist.scale());
 339   using boost::math::tools::max_value;
 340
 341   RealType min = 0.; // Minimum possible value is bottom of range of distribution.
 342   RealType max = max_value<RealType>();// Maximum possible value is top of range.
 343   // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
 344   // digits used to control how accurate to try to make the result.
 345   // To allow user to control accuracy versus speed,
 346   int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
 347   boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
 348   using boost::math::tools::newton_raphson_iterate;
 349   result =
 350     newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType>(dist, p), guess, min, max, get_digits, m);
 351    return result;
 352 } // quantile
 353
 354 template <class RealType, class Policy>
 355 inline RealType cdf(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)
 356 {
 357    BOOST_MATH_STD_USING  // for ADL of std functions.
 358
 359    RealType scale = c.dist.scale();
 360    RealType mean = c.dist.mean();
 361    RealType x = c.param;
 362    static const char* function = "boost::math::cdf(const complement(inverse_gaussian_distribution<%1%>&), %1%)";
 363    // infinite arguments not supported.
 364    //if((boost::math::isinf)(x))
 365    //{
 366    //  if(x < 0) return 1; // cdf complement -infinity is unity.
 367    //  return 0; // cdf complement +infinity is zero
 368    //}
 369    // These produce MSVC 4127 warnings, so the above used instead.
 370    //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
 371    //{ // cdf complement +infinity is zero.
 372    //  return 0;
 373    //}
 374    //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
 375    //{ // cdf complement -infinity is unity.
 376    //  return 1;
 377    //}
 378    RealType result = 0;
 379    if(false == detail::check_scale(function, scale, &result, Policy()))
 380       return result;
 381    if(false == detail::check_location(function, mean, &result, Policy()))
 382       return result;
 383    if(false == detail::check_x(function, x, &result, Policy()))
 384       return result;
 385
 386    normal_distribution<RealType> n01;
 387    RealType n0 = sqrt(scale / x);
 388    n0 *= ((x / mean) -1);
 389    RealType cdf_1 = cdf(complement(n01, n0));
 390
 391    RealType expfactor = exp(2 * scale / mean);
 392    RealType n3 = - sqrt(scale / x);
 393    n3 *= (x / mean) + 1;
 394
 395    //RealType n5 = +sqrt(scale/x) * ((x /mean) + 1); // note now positive sign.
 396    RealType n6 = cdf(complement(n01, +sqrt(scale/x) * ((x /mean) + 1)));
 397    // RealType n4 = cdf(n01, n3); // =
 398    result = cdf_1 - expfactor * n6;
 399    return result;
 400 } // cdf complement
 401
 402 template <class RealType, class Policy>
 403 inline RealType quantile(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)
 404 {
 405    BOOST_MATH_STD_USING  // for ADL of std functions
 406
 407    RealType scale = c.dist.scale();
 408    RealType mean = c.dist.mean();
 409    static const char* function = "boost::math::quantile(const complement(inverse_gaussian_distribution<%1%>&), %1%)";
 410    RealType result = 0;
 411    if(false == detail::check_scale(function, scale, &result, Policy()))
 412       return result;
 413    if(false == detail::check_location(function, mean, &result, Policy()))
 414       return result;
 415    RealType q = c.param;
 416    if(false == detail::check_probability(function, q, &result, Policy()))
 417       return result;
 418
 419    RealType guess = detail::guess_ig(q, mean, scale);
 420    // Complement.
 421    using boost::math::tools::max_value;
 422
 423   RealType min = 0.; // Minimum possible value is bottom of range of distribution.
 424   RealType max = max_value<RealType>();// Maximum possible value is top of range.
 425   // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
 426   // digits used to control how accurate to try to make the result.
 427   int get_digits = policies::digits<RealType, Policy>();
 428   boost::uintmax_t m = policies::get_max_root_iterations<Policy>();
 429   using boost::math::tools::newton_raphson_iterate;
 430   result =
 431     newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType>(c.dist, q), guess, min, max, get_digits, m);
 432    return result;
 433 } // quantile
 434
 435 template <class RealType, class Policy>
 436 inline RealType mean(const inverse_gaussian_distribution<RealType, Policy>& dist)
 437 { // aka mu
 438    return dist.mean();
 439 }
 440
 441 template <class RealType, class Policy>
 442 inline RealType scale(const inverse_gaussian_distribution<RealType, Policy>& dist)
 443 { // aka lambda
 444    return dist.scale();
 445 }
 446
 447 template <class RealType, class Policy>
 448 inline RealType shape(const inverse_gaussian_distribution<RealType, Policy>& dist)
 449 { // aka phi
 450    return dist.shape();
 451 }
 452
 453 template <class RealType, class Policy>
 454 inline RealType standard_deviation(const inverse_gaussian_distribution<RealType, Policy>& dist)
 455 {
 456   BOOST_MATH_STD_USING
 457   RealType scale = dist.scale();
 458   RealType mean = dist.mean();
 459   RealType result = sqrt(mean * mean * mean / scale);
 460   return result;
 461 }
 462
 463 template <class RealType, class Policy>
 464 inline RealType mode(const inverse_gaussian_distribution<RealType, Policy>& dist)
 465 {
 466   BOOST_MATH_STD_USING
 467   RealType scale = dist.scale();
 468   RealType  mean = dist.mean();
 469   RealType result = mean * (sqrt(1 + (9 * mean * mean)/(4 * scale * scale))
 470       - 3 * mean / (2 * scale));
 471   return result;
 472 }
 473
 474 template <class RealType, class Policy>
 475 inline RealType skewness(const inverse_gaussian_distribution<RealType, Policy>& dist)
 476 {
 477   BOOST_MATH_STD_USING
 478   RealType scale = dist.scale();
 479   RealType  mean = dist.mean();
 480   RealType result = 3 * sqrt(mean/scale);
 481   return result;
 482 }
 483
 484 template <class RealType, class Policy>
 485 inline RealType kurtosis(const inverse_gaussian_distribution<RealType, Policy>& dist)
 486 {
 487   RealType scale = dist.scale();
 488   RealType  mean = dist.mean();
 489   RealType result = 15 * mean / scale -3;
 490   return result;
 491 }
 492
 493 template <class RealType, class Policy>
 494 inline RealType kurtosis_excess(const inverse_gaussian_distribution<RealType, Policy>& dist)
 495 {
 496   RealType scale = dist.scale();
 497   RealType  mean = dist.mean();
 498   RealType result = 15 * mean / scale;
 499   return result;
 500 }
 501
 502 } // namespace math
 503 } // namespace boost
 504
 505 // This include must be at the end, *after* the accessors
 506 // for this distribution have been defined, in order to
 507 // keep compilers that support two-phase lookup happy.
 508 #include <boost/math/distributions/detail/derived_accessors.hpp>
 509
 510 #endif // BOOST_STATS_INVERSE_GAUSSIAN_HPP
 511
 512