Chris@16: //  Copyright John Maddock 2010.
Chris@16: //  Copyright Paul A. Bristow 2010.
Chris@16: 
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_STATS_INVERSE_GAUSSIAN_HPP
Chris@16: #define BOOST_STATS_INVERSE_GAUSSIAN_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma warning(disable: 4512) // assignment operator could not be generated
Chris@16: #endif
Chris@16: 
Chris@16: // http://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution
Chris@16: // http://mathworld.wolfram.com/InverseGaussianDistribution.html
Chris@16: 
Chris@16: // The normal-inverse Gaussian distribution
Chris@16: // also called the Wald distribution (some sources limit this to when mean = 1).
Chris@16: 
Chris@16: // It is the continuous probability distribution
Chris@16: // that is defined as the normal variance-mean mixture where the mixing density is the 
Chris@16: // inverse Gaussian distribution. The tails of the distribution decrease more slowly
Chris@16: // than the normal distribution. It is therefore suitable to model phenomena
Chris@16: // where numerically large values are more probable than is the case for the normal distribution.
Chris@16: 
Chris@16: // The Inverse Gaussian distribution was first studied in relationship to Brownian motion.
Chris@16: // In 1956 M.C.K. Tweedie used the name 'Inverse Gaussian' because there is an inverse 
Chris@16: // relationship between the time to cover a unit distance and distance covered in unit time.
Chris@16: 
Chris@16: // Examples are returns from financial assets and turbulent wind speeds. 
Chris@16: // The normal-inverse Gaussian distributions form
Chris@16: // a subclass of the generalised hyperbolic distributions.
Chris@16: 
Chris@16: // See also
Chris@16: 
Chris@16: // http://en.wikipedia.org/wiki/Normal_distribution
Chris@16: // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
Chris@16: // Also:
Chris@16: // Weisstein, Eric W. "Normal Distribution."
Chris@16: // From MathWorld--A Wolfram Web Resource.
Chris@16: // http://mathworld.wolfram.com/NormalDistribution.html
Chris@16: 
Chris@16: // http://www.jstatsoft.org/v26/i04/paper General class of inverse Gaussian distributions.
Chris@16: // ig package - withdrawn but at http://cran.r-project.org/src/contrib/Archive/ig/
Chris@16: 
Chris@16: // http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/SuppDists/html/inverse_gaussian.html
Chris@16: // R package for dinverse_gaussian, ...
Chris@16: 
Chris@16: // http://www.statsci.org/s/inverse_gaussian.s  and http://www.statsci.org/s/inverse_gaussian.html
Chris@16: 
Chris@16: //#include <boost/math/distributions/fwd.hpp>
Chris@16: #include <boost/math/special_functions/erf.hpp> // for erf/erfc.
Chris@16: #include <boost/math/distributions/complement.hpp>
Chris@16: #include <boost/math/distributions/detail/common_error_handling.hpp>
Chris@16: #include <boost/math/distributions/normal.hpp>
Chris@16: #include <boost/math/distributions/gamma.hpp> // for gamma function
Chris@16: // using boost::math::gamma_p;
Chris@16: 
Chris@16: #include <boost/math/tools/tuple.hpp>
Chris@16: //using std::tr1::tuple;
Chris@16: //using std::tr1::make_tuple;
Chris@16: #include <boost/math/tools/roots.hpp>
Chris@16: //using boost::math::tools::newton_raphson_iterate;
Chris@16: 
Chris@16: #include <utility>
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: template <class RealType = double, class Policy = policies::policy<> >
Chris@16: class inverse_gaussian_distribution
Chris@16: {
Chris@16: public:
Chris@16:    typedef RealType value_type;
Chris@16:    typedef Policy policy_type;
Chris@16: 
Chris@16:    inverse_gaussian_distribution(RealType l_mean = 1, RealType l_scale = 1)
Chris@16:       : m_mean(l_mean), m_scale(l_scale)
Chris@16:    { // Default is a 1,1 inverse_gaussian distribution.
Chris@16:      static const char* function = "boost::math::inverse_gaussian_distribution<%1%>::inverse_gaussian_distribution";
Chris@16: 
Chris@16:      RealType result;
Chris@16:      detail::check_scale(function, l_scale, &result, Policy());
Chris@16:      detail::check_location(function, l_mean, &result, Policy());
Chris@16:    }
Chris@16: 
Chris@16:    RealType mean()const
Chris@16:    { // alias for location.
Chris@16:       return m_mean; // aka mu
Chris@16:    }
Chris@16: 
Chris@16:    // Synonyms, provided to allow generic use of find_location and find_scale.
Chris@16:    RealType location()const
Chris@16:    { // location, aka mu.
Chris@16:       return m_mean;
Chris@16:    }
Chris@16:    RealType scale()const
Chris@16:    { // scale, aka lambda.
Chris@16:       return m_scale;
Chris@16:    }
Chris@16: 
Chris@16:    RealType shape()const
Chris@16:    { // shape, aka phi = lambda/mu.
Chris@16:       return m_scale / m_mean;
Chris@16:    }
Chris@16: 
Chris@16: private:
Chris@16:    //
Chris@16:    // Data members:
Chris@16:    //
Chris@16:    RealType m_mean;  // distribution mean or location, aka mu.
Chris@16:    RealType m_scale;    // distribution standard deviation or scale, aka lambda.
Chris@16: }; // class normal_distribution
Chris@16: 
Chris@16: typedef inverse_gaussian_distribution<double> inverse_gaussian;
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline const std::pair<RealType, RealType> range(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)
Chris@16: { // Range of permissible values for random variable x, zero to max.
Chris@16:    using boost::math::tools::max_value;
Chris@16:    return std::pair<RealType, RealType>(static_cast<RealType>(0.), max_value<RealType>()); // - to + max value.
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline const std::pair<RealType, RealType> support(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)
Chris@16: { // Range of supported values for random variable x, zero to max.
Chris@16:   // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
Chris@16:    using boost::math::tools::max_value;
Chris@16:    return std::pair<RealType, RealType>(static_cast<RealType>(0.),  max_value<RealType>()); // - to + max value.
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType pdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
Chris@16: { // Probability Density Function
Chris@16:    BOOST_MATH_STD_USING  // for ADL of std functions
Chris@16: 
Chris@16:    RealType scale = dist.scale();
Chris@16:    RealType mean = dist.mean();
Chris@16:    RealType result = 0;
Chris@16:    static const char* function = "boost::math::pdf(const inverse_gaussian_distribution<%1%>&, %1%)";
Chris@16:    if(false == detail::check_scale(function, scale, &result, Policy()))
Chris@16:    {
Chris@16:       return result;
Chris@16:    }
Chris@16:    if(false == detail::check_location(function, mean, &result, Policy()))
Chris@16:    {
Chris@16:       return result;
Chris@16:    }
Chris@16:    if(false == detail::check_positive_x(function, x, &result, Policy()))
Chris@16:    {
Chris@16:       return result;
Chris@16:    }
Chris@16: 
Chris@16:    if (x == 0)
Chris@16:    {
Chris@16:      return 0; // Convenient, even if not defined mathematically.
Chris@16:    }
Chris@16: 
Chris@16:    result =
Chris@16:      sqrt(scale / (constants::two_pi<RealType>() * x * x * x))
Chris@16:     * exp(-scale * (x - mean) * (x - mean) / (2 * x * mean * mean));
Chris@16:    return result;
Chris@16: } // pdf
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType cdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
Chris@16: { // Cumulative Density Function.
Chris@16:    BOOST_MATH_STD_USING  // for ADL of std functions.
Chris@16: 
Chris@16:    RealType scale = dist.scale();
Chris@16:    RealType mean = dist.mean();
Chris@16:    static const char* function = "boost::math::cdf(const inverse_gaussian_distribution<%1%>&, %1%)";
Chris@16:    RealType result = 0;
Chris@16:    if(false == detail::check_scale(function, scale, &result, Policy()))
Chris@16:    {
Chris@16:       return result;
Chris@16:    }
Chris@16:    if(false == detail::check_location(function, mean, &result, Policy()))
Chris@16:    {
Chris@16:       return result;
Chris@16:    }
Chris@16:    if(false == detail::check_positive_x(function, x, &result, Policy()))
Chris@16:    {
Chris@16:      return result;
Chris@16:    }
Chris@16:    if (x == 0)
Chris@16:    {
Chris@16:      return 0; // Convenient, even if not defined mathematically.
Chris@16:    }
Chris@16:    // Problem with this formula for large scale > 1000 or small x, 
Chris@16:    //result = 0.5 * (erf(sqrt(scale / x) * ((x / mean) - 1) / constants::root_two<RealType>(), Policy()) + 1)
Chris@16:    //  + exp(2 * scale / mean) / 2 
Chris@16:    //  * (1 - erf(sqrt(scale / x) * (x / mean + 1) / constants::root_two<RealType>(), Policy()));
Chris@16:    // so use normal distribution version:
Chris@16:    // Wikipedia CDF equation http://en.wikipedia.org/wiki/Inverse_Gaussian_distribution.
Chris@16: 
Chris@16:    normal_distribution<RealType> n01;
Chris@16: 
Chris@16:    RealType n0 = sqrt(scale / x);
Chris@16:    n0 *= ((x / mean) -1);
Chris@16:    RealType n1 = cdf(n01, n0);
Chris@16:    RealType expfactor = exp(2 * scale / mean);
Chris@16:    RealType n3 = - sqrt(scale / x);
Chris@16:    n3 *= (x / mean) + 1;
Chris@16:    RealType n4 = cdf(n01, n3);
Chris@16:    result = n1 + expfactor * n4;
Chris@16:    return result;
Chris@16: } // cdf
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: struct inverse_gaussian_quantile_functor
Chris@16: { 
Chris@16: 
Chris@16:   inverse_gaussian_quantile_functor(const boost::math::inverse_gaussian_distribution<RealType, Policy> dist, RealType const& p)
Chris@16:     : distribution(dist), prob(p)
Chris@16:   {
Chris@16:   }
Chris@16:   boost::math::tuple<RealType, RealType> operator()(RealType const& x)
Chris@16:   {
Chris@16:     RealType c = cdf(distribution, x);
Chris@16:     RealType fx = c - prob;  // Difference cdf - value - to minimize.
Chris@16:     RealType dx = pdf(distribution, x); // pdf is 1st derivative.
Chris@16:     // return both function evaluation difference f(x) and 1st derivative f'(x).
Chris@16:     return boost::math::make_tuple(fx, dx);
Chris@16:   }
Chris@16:   private:
Chris@16:   const boost::math::inverse_gaussian_distribution<RealType, Policy> distribution;
Chris@16:   RealType prob; 
Chris@16: };
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: struct inverse_gaussian_quantile_complement_functor
Chris@16: { 
Chris@16:     inverse_gaussian_quantile_complement_functor(const boost::math::inverse_gaussian_distribution<RealType, Policy> dist, RealType const& p)
Chris@16:     : distribution(dist), prob(p)
Chris@16:   {
Chris@16:   }
Chris@16:   boost::math::tuple<RealType, RealType> operator()(RealType const& x)
Chris@16:   {
Chris@16:     RealType c = cdf(complement(distribution, x));
Chris@16:     RealType fx = c - prob;  // Difference cdf - value - to minimize.
Chris@16:     RealType dx = -pdf(distribution, x); // pdf is 1st derivative.
Chris@16:     // return both function evaluation difference f(x) and 1st derivative f'(x).
Chris@16:     //return std::tr1::make_tuple(fx, dx); if available.
Chris@16:     return boost::math::make_tuple(fx, dx);
Chris@16:   }
Chris@16:   private:
Chris@16:   const boost::math::inverse_gaussian_distribution<RealType, Policy> distribution;
Chris@16:   RealType prob; 
Chris@16: };
Chris@16: 
Chris@16: namespace detail
Chris@16: {
Chris@16:   template <class RealType>
Chris@16:   inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1)
Chris@16:   { // guess at random variate value x for inverse gaussian quantile.
Chris@16:       BOOST_MATH_STD_USING
Chris@16:       using boost::math::policies::policy;
Chris@16:       // Error type.
Chris@16:       using boost::math::policies::overflow_error;
Chris@16:       // Action.
Chris@16:       using boost::math::policies::ignore_error;
Chris@16: 
Chris@16:       typedef policy<
Chris@16:         overflow_error<ignore_error> // Ignore overflow (return infinity)
Chris@16:       > no_overthrow_policy;
Chris@16: 
Chris@16:     RealType x; // result is guess at random variate value x.
Chris@16:     RealType phi = lambda / mu;
Chris@16:     if (phi > 2.)
Chris@16:     { // Big phi, so starting to look like normal Gaussian distribution.
Chris@16:       //    x=(qnorm(p,0,1,true,false) - 0.5 * sqrt(mu/lambda)) / sqrt(lambda/mu);
Chris@16:       // Whitmore, G.A. and Yalovsky, M.
Chris@16:       // A normalising logarithmic transformation for inverse Gaussian random variables,
Chris@16:       // Technometrics 20-2, 207-208 (1978), but using expression from
Chris@16:       // V Seshadri, Inverse Gaussian distribution (1998) ISBN 0387 98618 9, page 6.
Chris@16:  
Chris@16:       normal_distribution<RealType, no_overthrow_policy> n01;
Chris@16:       x = mu * exp(quantile(n01, p) / sqrt(phi) - 1/(2 * phi));
Chris@16:      }
Chris@16:     else
Chris@16:     { // phi < 2 so much less symmetrical with long tail,
Chris@16:       // so use gamma distribution as an approximation.
Chris@16:       using boost::math::gamma_distribution;
Chris@16: 
Chris@16:       // Define the distribution, using gamma_nooverflow:
Chris@16:       typedef gamma_distribution<RealType, no_overthrow_policy> gamma_nooverflow;
Chris@16: 
Chris@16:       gamma_nooverflow g(static_cast<RealType>(0.5), static_cast<RealType>(1.));
Chris@16: 
Chris@16:       // gamma_nooverflow g(static_cast<RealType>(0.5), static_cast<RealType>(1.));
Chris@16:       // R qgamma(0.2, 0.5, 1)  0.0320923
Chris@16:       RealType qg = quantile(complement(g, p));
Chris@16:       //RealType qg1 = qgamma(1.- p, 0.5, 1.0, true, false);
Chris@16:       x = lambda / (qg * 2);
Chris@16:       // 
Chris@16:       if (x > mu/2) // x > mu /2?
Chris@16:       { // x too large for the gamma approximation to work well.
Chris@16:         //x = qgamma(p, 0.5, 1.0); // qgamma(0.270614, 0.5, 1) = 0.05983807
Chris@16:         RealType q = quantile(g, p);
Chris@16:        // x = mu * exp(q * static_cast<RealType>(0.1));  // Said to improve at high p
Chris@16:        // x = mu * x;  // Improves at high p?
Chris@16:         x = mu * exp(q / sqrt(phi) - 1/(2 * phi));
Chris@16:       }
Chris@16:     }
Chris@16:     return x;
Chris@16:   }  // guess_ig
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType quantile(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& p)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // for ADL of std functions.
Chris@16:    // No closed form exists so guess and use Newton Raphson iteration.
Chris@16: 
Chris@16:    RealType mean = dist.mean();
Chris@16:    RealType scale = dist.scale();
Chris@16:    static const char* function = "boost::math::quantile(const inverse_gaussian_distribution<%1%>&, %1%)";
Chris@16: 
Chris@16:    RealType result = 0;
Chris@16:    if(false == detail::check_scale(function, scale, &result, Policy()))
Chris@16:       return result;
Chris@16:    if(false == detail::check_location(function, mean, &result, Policy()))
Chris@16:       return result;
Chris@16:    if(false == detail::check_probability(function, p, &result, Policy()))
Chris@16:       return result;
Chris@16:    if (p == 0)
Chris@16:    {
Chris@16:      return 0; // Convenient, even if not defined mathematically?
Chris@16:    }
Chris@16:    if (p == 1)
Chris@16:    { // overflow 
Chris@16:       result = policies::raise_overflow_error<RealType>(function,
Chris@16:         "probability parameter is 1, but must be < 1!", Policy());
Chris@16:       return result; // std::numeric_limits<RealType>::infinity();
Chris@16:    }
Chris@16: 
Chris@16:   RealType guess = detail::guess_ig(p, dist.mean(), dist.scale());
Chris@16:   using boost::math::tools::max_value;
Chris@16: 
Chris@16:   RealType min = 0.; // Minimum possible value is bottom of range of distribution.
Chris@16:   RealType max = max_value<RealType>();// Maximum possible value is top of range. 
Chris@16:   // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
Chris@16:   // digits used to control how accurate to try to make the result.
Chris@16:   // To allow user to control accuracy versus speed,
Chris@16:   int get_digits = policies::digits<RealType, Policy>();// get digits from policy, 
Chris@16:   boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
Chris@16:   using boost::math::tools::newton_raphson_iterate;
Chris@16:   result =
Chris@16:     newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType, Policy>(dist, p), guess, min, max, get_digits, m);
Chris@16:    return result;
Chris@16: } // quantile
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType cdf(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // for ADL of std functions.
Chris@16: 
Chris@16:    RealType scale = c.dist.scale();
Chris@16:    RealType mean = c.dist.mean();
Chris@16:    RealType x = c.param;
Chris@16:    static const char* function = "boost::math::cdf(const complement(inverse_gaussian_distribution<%1%>&), %1%)";
Chris@16:    // infinite arguments not supported.
Chris@16:    //if((boost::math::isinf)(x))
Chris@16:    //{
Chris@16:    //  if(x < 0) return 1; // cdf complement -infinity is unity.
Chris@16:    //  return 0; // cdf complement +infinity is zero
Chris@16:    //}
Chris@16:    // These produce MSVC 4127 warnings, so the above used instead.
Chris@16:    //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
Chris@16:    //{ // cdf complement +infinity is zero.
Chris@16:    //  return 0;
Chris@16:    //}
Chris@16:    //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
Chris@16:    //{ // cdf complement -infinity is unity.
Chris@16:    //  return 1;
Chris@16:    //}
Chris@16:    RealType result = 0;
Chris@16:    if(false == detail::check_scale(function, scale, &result, Policy()))
Chris@16:       return result;
Chris@16:    if(false == detail::check_location(function, mean, &result, Policy()))
Chris@16:       return result;
Chris@16:    if(false == detail::check_positive_x(function, x, &result, Policy()))
Chris@16:       return result;
Chris@16: 
Chris@16:    normal_distribution<RealType> n01;
Chris@16:    RealType n0 = sqrt(scale / x);
Chris@16:    n0 *= ((x / mean) -1);
Chris@16:    RealType cdf_1 = cdf(complement(n01, n0));
Chris@16: 
Chris@16:    RealType expfactor = exp(2 * scale / mean);
Chris@16:    RealType n3 = - sqrt(scale / x);
Chris@16:    n3 *= (x / mean) + 1;
Chris@16: 
Chris@16:    //RealType n5 = +sqrt(scale/x) * ((x /mean) + 1); // note now positive sign.
Chris@16:    RealType n6 = cdf(complement(n01, +sqrt(scale/x) * ((x /mean) + 1)));
Chris@16:    // RealType n4 = cdf(n01, n3); // = 
Chris@16:    result = cdf_1 - expfactor * n6; 
Chris@16:    return result;
Chris@16: } // cdf complement
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType quantile(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // for ADL of std functions
Chris@16: 
Chris@16:    RealType scale = c.dist.scale();
Chris@16:    RealType mean = c.dist.mean();
Chris@16:    static const char* function = "boost::math::quantile(const complement(inverse_gaussian_distribution<%1%>&), %1%)";
Chris@16:    RealType result = 0;
Chris@16:    if(false == detail::check_scale(function, scale, &result, Policy()))
Chris@16:       return result;
Chris@16:    if(false == detail::check_location(function, mean, &result, Policy()))
Chris@16:       return result;
Chris@16:    RealType q = c.param;
Chris@16:    if(false == detail::check_probability(function, q, &result, Policy()))
Chris@16:       return result;
Chris@16: 
Chris@16:    RealType guess = detail::guess_ig(q, mean, scale);
Chris@16:    // Complement.
Chris@16:    using boost::math::tools::max_value;
Chris@16: 
Chris@16:   RealType min = 0.; // Minimum possible value is bottom of range of distribution.
Chris@16:   RealType max = max_value<RealType>();// Maximum possible value is top of range. 
Chris@16:   // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
Chris@16:   // digits used to control how accurate to try to make the result.
Chris@16:   int get_digits = policies::digits<RealType, Policy>();
Chris@16:   boost::uintmax_t m = policies::get_max_root_iterations<Policy>();
Chris@16:   using boost::math::tools::newton_raphson_iterate;
Chris@16:   result =
Chris@16:     newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType, Policy>(c.dist, q), guess, min, max, get_digits, m);
Chris@16:    return result;
Chris@16: } // quantile
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType mean(const inverse_gaussian_distribution<RealType, Policy>& dist)
Chris@16: { // aka mu
Chris@16:    return dist.mean();
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType scale(const inverse_gaussian_distribution<RealType, Policy>& dist)
Chris@16: { // aka lambda
Chris@16:    return dist.scale();
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType shape(const inverse_gaussian_distribution<RealType, Policy>& dist)
Chris@16: { // aka phi
Chris@16:    return dist.shape();
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType standard_deviation(const inverse_gaussian_distribution<RealType, Policy>& dist)
Chris@16: {
Chris@16:   BOOST_MATH_STD_USING
Chris@16:   RealType scale = dist.scale();
Chris@16:   RealType mean = dist.mean();
Chris@16:   RealType result = sqrt(mean * mean * mean / scale);
Chris@16:   return result;
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType mode(const inverse_gaussian_distribution<RealType, Policy>& dist)
Chris@16: {
Chris@16:   BOOST_MATH_STD_USING
Chris@16:   RealType scale = dist.scale();
Chris@16:   RealType  mean = dist.mean();
Chris@16:   RealType result = mean * (sqrt(1 + (9 * mean * mean)/(4 * scale * scale)) 
Chris@16:       - 3 * mean / (2 * scale));
Chris@16:   return result;
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType skewness(const inverse_gaussian_distribution<RealType, Policy>& dist)
Chris@16: {
Chris@16:   BOOST_MATH_STD_USING
Chris@16:   RealType scale = dist.scale();
Chris@16:   RealType  mean = dist.mean();
Chris@16:   RealType result = 3 * sqrt(mean/scale);
Chris@16:   return result;
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType kurtosis(const inverse_gaussian_distribution<RealType, Policy>& dist)
Chris@16: {
Chris@16:   RealType scale = dist.scale();
Chris@16:   RealType  mean = dist.mean();
Chris@16:   RealType result = 15 * mean / scale -3;
Chris@16:   return result;
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType kurtosis_excess(const inverse_gaussian_distribution<RealType, Policy>& dist)
Chris@16: {
Chris@16:   RealType scale = dist.scale();
Chris@16:   RealType  mean = dist.mean();
Chris@16:   RealType result = 15 * mean / scale;
Chris@16:   return result;
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: // This include must be at the end, *after* the accessors
Chris@16: // for this distribution have been defined, in order to
Chris@16: // keep compilers that support two-phase lookup happy.
Chris@16: #include <boost/math/distributions/detail/derived_accessors.hpp>
Chris@16: 
Chris@16: #endif // BOOST_STATS_INVERSE_GAUSSIAN_HPP
Chris@16: 
Chris@16: