SV platform dependency builds

//  Copyright John Maddock 2010.

//  Copyright Paul A. Bristow 2010.

//  Use, modification and distribution are subject to the

//  Boost Software License, Version 1.0. (See accompanying file

//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_STATS_INVERSE_GAUSSIAN_HPP

#define BOOST_STATS_INVERSE_GAUSSIAN_HPP

#ifdef _MSC_VER

#pragma warning(disable: 4512) // assignment operator could not be generated

#endif

// http://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution

// http://mathworld.wolfram.com/InverseGaussianDistribution.html

// The normal-inverse Gaussian distribution

// also called the Wald distribution (some sources limit this to when mean = 1).

// It is the continuous probability distribution

// that is defined as the normal variance-mean mixture where the mixing density is the 

// inverse Gaussian distribution. The tails of the distribution decrease more slowly

// than the normal distribution. It is therefore suitable to model phenomena

// where numerically large values are more probable than is the case for the normal distribution.

// The Inverse Gaussian distribution was first studied in relationship to Brownian motion.

// In 1956 M.C.K. Tweedie used the name 'Inverse Gaussian' because there is an inverse 

// relationship between the time to cover a unit distance and distance covered in unit time.

// Examples are returns from financial assets and turbulent wind speeds. 

// The normal-inverse Gaussian distributions form

// a subclass of the generalised hyperbolic distributions.

// See also

// http://en.wikipedia.org/wiki/Normal_distribution

// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm

// Also:

// Weisstein, Eric W. "Normal Distribution."

// From MathWorld--A Wolfram Web Resource.

// http://mathworld.wolfram.com/NormalDistribution.html

// http://www.jstatsoft.org/v26/i04/paper General class of inverse Gaussian distributions.

// ig package - withdrawn but at http://cran.r-project.org/src/contrib/Archive/ig/

// http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/SuppDists/html/inverse_gaussian.html

// R package for dinverse_gaussian, ...

// http://www.statsci.org/s/inverse_gaussian.s  and http://www.statsci.org/s/inverse_gaussian.html

//#include <boost/math/distributions/fwd.hpp>

#include <boost/math/special_functions/erf.hpp> // for erf/erfc.

#include <boost/math/distributions/complement.hpp>

#include <boost/math/distributions/detail/common_error_handling.hpp>

#include <boost/math/distributions/normal.hpp>

#include <boost/math/distributions/gamma.hpp> // for gamma function

// using boost::math::gamma_p;

#include <boost/math/tools/tuple.hpp>

//using std::tr1::tuple;

//using std::tr1::make_tuple;

#include <boost/math/tools/roots.hpp>

//using boost::math::tools::newton_raphson_iterate;

#include <utility>

namespace boost{ namespace math{

template <class RealType = double, class Policy = policies::policy<> >

class inverse_gaussian_distribution

{

public:

   typedef RealType value_type;

   typedef Policy policy_type;

   inverse_gaussian_distribution(RealType l_mean = 1, RealType l_scale = 1)

      : m_mean(l_mean), m_scale(l_scale)

   { // Default is a 1,1 inverse_gaussian distribution.

     static const char* function = "boost::math::inverse_gaussian_distribution<%1%>::inverse_gaussian_distribution";

     RealType result;

     detail::check_scale(function, l_scale, &result, Policy());

     detail::check_location(function, l_mean, &result, Policy());

     detail::check_x_gt0(function, l_mean, &result, Policy());

   }

   RealType mean()const

   { // alias for location.

      return m_mean; // aka mu

   }

   // Synonyms, provided to allow generic use of find_location and find_scale.

   RealType location()const

   { // location, aka mu.

      return m_mean;

   }

   RealType scale()const

   { // scale, aka lambda.

      return m_scale;

   }

   RealType shape()const

   { // shape, aka phi = lambda/mu.

      return m_scale / m_mean;

   }

private:

   //

   // Data members:

   //

   RealType m_mean;  // distribution mean or location, aka mu.

   RealType m_scale;    // distribution standard deviation or scale, aka lambda.

}; // class normal_distribution

typedef inverse_gaussian_distribution<double> inverse_gaussian;

template <class RealType, class Policy>

inline const std::pair<RealType, RealType> range(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)

{ // Range of permissible values for random variable x, zero to max.

   using boost::math::tools::max_value;

   return std::pair<RealType, RealType>(static_cast<RealType>(0.), max_value<RealType>()); // - to + max value.

}

template <class RealType, class Policy>

inline const std::pair<RealType, RealType> support(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)

{ // Range of supported values for random variable x, zero to max.

  // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.

   using boost::math::tools::max_value;

   return std::pair<RealType, RealType>(static_cast<RealType>(0.),  max_value<RealType>()); // - to + max value.

}

template <class RealType, class Policy>

inline RealType pdf(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::pdf(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.

   }

   result =

     sqrt(scale / (constants::two_pi<RealType>() * x * x * x))

    * exp(-scale * (x - mean) * (x - mean) / (2 * x * mean * mean));

   return result;

} // pdf

template <class RealType, class Policy>

inline RealType cdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)

{ // Cumulative Density Function.

   BOOST_MATH_STD_USING  // for ADL of std functions.

   RealType scale = dist.scale();

   RealType mean = dist.mean();

   static const char* function = "boost::math::cdf(const inverse_gaussian_distribution<%1%>&, %1%)";

   RealType result = 0;

   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.

   }

   // Problem with this formula for large scale > 1000 or small x, 

   //result = 0.5 * (erf(sqrt(scale / x) * ((x / mean) - 1) / constants::root_two<RealType>(), Policy()) + 1)

   //  + exp(2 * scale / mean) / 2 

   //  * (1 - erf(sqrt(scale / x) * (x / mean + 1) / constants::root_two<RealType>(), Policy()));

   // so use normal distribution version:

   // Wikipedia CDF equation http://en.wikipedia.org/wiki/Inverse_Gaussian_distribution.

   normal_distribution<RealType> n01;

   RealType n0 = sqrt(scale / x);

   n0 *= ((x / mean) -1);

   RealType n1 = cdf(n01, n0);

   RealType expfactor = exp(2 * scale / mean);

   RealType n3 = - sqrt(scale / x);

   n3 *= (x / mean) + 1;

   RealType n4 = cdf(n01, n3);

   result = n1 + expfactor * n4;

   return result;

} // cdf

template <class RealType, class Policy>

struct inverse_gaussian_quantile_functor

{ 

  inverse_gaussian_quantile_functor(const boost::math::inverse_gaussian_distribution<RealType, Policy> dist, RealType const& p)

    : distribution(dist), prob(p)

  {

  }

  boost::math::tuple<RealType, RealType> operator()(RealType const& x)

  {

    RealType c = cdf(distribution, x);

    RealType fx = c - prob;  // Difference cdf - value - to minimize.

    RealType dx = pdf(distribution, x); // pdf is 1st derivative.

    // return both function evaluation difference f(x) and 1st derivative f'(x).

    return boost::math::make_tuple(fx, dx);

  }

  private:

  const boost::math::inverse_gaussian_distribution<RealType, Policy> distribution;

  RealType prob; 

};

template <class RealType, class Policy>

struct inverse_gaussian_quantile_complement_functor

{ 

    inverse_gaussian_quantile_complement_functor(const boost::math::inverse_gaussian_distribution<RealType, Policy> dist, RealType const& p)

    : distribution(dist), prob(p)

  {

  }

  boost::math::tuple<RealType, RealType> operator()(RealType const& x)

  {

    RealType c = cdf(complement(distribution, x));

    RealType fx = c - prob;  // Difference cdf - value - to minimize.

    RealType dx = -pdf(distribution, x); // pdf is 1st derivative.

    // return both function evaluation difference f(x) and 1st derivative f'(x).

    //return std::tr1::make_tuple(fx, dx); if available.

    return boost::math::make_tuple(fx, dx);

  }

  private:

  const boost::math::inverse_gaussian_distribution<RealType, Policy> distribution;

  RealType prob; 

};

namespace detail

{

  template <class RealType>

  inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1)

  { // guess at random variate value x for inverse gaussian quantile.

      BOOST_MATH_STD_USING

      using boost::math::policies::policy;

      // Error type.

      using boost::math::policies::overflow_error;

      // Action.

      using boost::math::policies::ignore_error;

      typedef policy<

        overflow_error<ignore_error> // Ignore overflow (return infinity)

      > no_overthrow_policy;

    RealType x; // result is guess at random variate value x.

    RealType phi = lambda / mu;

    if (phi > 2.)

    { // Big phi, so starting to look like normal Gaussian distribution.

      //    x=(qnorm(p,0,1,true,false) - 0.5 * sqrt(mu/lambda)) / sqrt(lambda/mu);

      // Whitmore, G.A. and Yalovsky, M.

      // A normalising logarithmic transformation for inverse Gaussian random variables,

      // Technometrics 20-2, 207-208 (1978), but using expression from

      // V Seshadri, Inverse Gaussian distribution (1998) ISBN 0387 98618 9, page 6.

      normal_distribution<RealType, no_overthrow_policy> n01;

      x = mu * exp(quantile(n01, p) / sqrt(phi) - 1/(2 * phi));

     }

    else

    { // phi < 2 so much less symmetrical with long tail,

      // so use gamma distribution as an approximation.

      using boost::math::gamma_distribution;

      // Define the distribution, using gamma_nooverflow:

      typedef gamma_distribution<RealType, no_overthrow_policy> gamma_nooverflow;

      gamma_nooverflow g(static_cast<RealType>(0.5), static_cast<RealType>(1.));

      // gamma_nooverflow g(static_cast<RealType>(0.5), static_cast<RealType>(1.));

      // R qgamma(0.2, 0.5, 1)  0.0320923

      RealType qg = quantile(complement(g, p));

      //RealType qg1 = qgamma(1.- p, 0.5, 1.0, true, false);

      x = lambda / (qg * 2);

      // 

      if (x > mu/2) // x > mu /2?

      { // x too large for the gamma approximation to work well.

        //x = qgamma(p, 0.5, 1.0); // qgamma(0.270614, 0.5, 1) = 0.05983807

        RealType q = quantile(g, p);

       // x = mu * exp(q * static_cast<RealType>(0.1));  // Said to improve at high p

       // x = mu * x;  // Improves at high p?

        x = mu * exp(q / sqrt(phi) - 1/(2 * phi));

      }

    }

    return x;

  }  // guess_ig

} // namespace detail

template <class RealType, class Policy>

inline RealType quantile(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& p)

{

   BOOST_MATH_STD_USING  // for ADL of std functions.

   // No closed form exists so guess and use Newton Raphson iteration.

   RealType mean = dist.mean();

   RealType scale = dist.scale();

   static const char* function = "boost::math::quantile(const inverse_gaussian_distribution<%1%>&, %1%)";

   RealType result = 0;

   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_probability(function, p, &result, Policy()))

      return result;

   if (p == 0)

   {

     return 0; // Convenient, even if not defined mathematically?

   }

   if (p == 1)

   { // overflow 

      result = policies::raise_overflow_error<RealType>(function,

        "probability parameter is 1, but must be < 1!", Policy());

      return result; // std::numeric_limits<RealType>::infinity();

   }

  RealType guess = detail::guess_ig(p, dist.mean(), dist.scale());

  using boost::math::tools::max_value;

  RealType min = 0.; // Minimum possible value is bottom of range of distribution.

  RealType max = max_value<RealType>();// Maximum possible value is top of range. 

  // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.

  // digits used to control how accurate to try to make the result.

  // To allow user to control accuracy versus speed,

  int get_digits = policies::digits<RealType, Policy>();// get digits from policy, 

  boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.

  using boost::math::tools::newton_raphson_iterate;

  result =

    newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType, Policy>(dist, p), guess, min, max, get_digits, m);

   return result;

} // quantile

template <class RealType, class Policy>

inline RealType cdf(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)

{

   BOOST_MATH_STD_USING  // for ADL of std functions.

   RealType scale = c.dist.scale();

   RealType mean = c.dist.mean();

   RealType x = c.param;

   static const char* function = "boost::math::cdf(const complement(inverse_gaussian_distribution<%1%>&), %1%)";

   // infinite arguments not supported.

   //if((boost::math::isinf)(x))

   //{

   //  if(x < 0) return 1; // cdf complement -infinity is unity.

   //  return 0; // cdf complement +infinity is zero

   //}

   // These produce MSVC 4127 warnings, so the above used instead.

   //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())

   //{ // cdf complement +infinity is zero.

   //  return 0;

   //}

   //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())

   //{ // cdf complement -infinity is unity.

   //  return 1;

   //}

   RealType result = 0;

   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;

   normal_distribution<RealType> n01;

   RealType n0 = sqrt(scale / x);

   n0 *= ((x / mean) -1);

   RealType cdf_1 = cdf(complement(n01, n0));

   RealType expfactor = exp(2 * scale / mean);

   RealType n3 = - sqrt(scale / x);

   n3 *= (x / mean) + 1;

   //RealType n5 = +sqrt(scale/x) * ((x /mean) + 1); // note now positive sign.

   RealType n6 = cdf(complement(n01, +sqrt(scale/x) * ((x /mean) + 1)));

   // RealType n4 = cdf(n01, n3); // = 

   result = cdf_1 - expfactor * n6; 

   return result;

} // cdf complement

template <class RealType, class Policy>

inline RealType quantile(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)

{

   BOOST_MATH_STD_USING  // for ADL of std functions

   RealType scale = c.dist.scale();

   RealType mean = c.dist.mean();

   static const char* function = "boost::math::quantile(const complement(inverse_gaussian_distribution<%1%>&), %1%)";

   RealType result = 0;

   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;

   RealType q = c.param;

   if(false == detail::check_probability(function, q, &result, Policy()))

      return result;

   RealType guess = detail::guess_ig(q, mean, scale);

   // Complement.

   using boost::math::tools::max_value;

  RealType min = 0.; // Minimum possible value is bottom of range of distribution.

  RealType max = max_value<RealType>();// Maximum possible value is top of range. 

  // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.

  // digits used to control how accurate to try to make the result.

  int get_digits = policies::digits<RealType, Policy>();

  boost::uintmax_t m = policies::get_max_root_iterations<Policy>();

  using boost::math::tools::newton_raphson_iterate;

  result =

    newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType, Policy>(c.dist, q), guess, min, max, get_digits, m);

   return result;

} // quantile

template <class RealType, class Policy>

inline RealType mean(const inverse_gaussian_distribution<RealType, Policy>& dist)

{ // aka mu

   return dist.mean();

}

template <class RealType, class Policy>

inline RealType scale(const inverse_gaussian_distribution<RealType, Policy>& dist)

{ // aka lambda

   return dist.scale();

}

template <class RealType, class Policy>

inline RealType shape(const inverse_gaussian_distribution<RealType, Policy>& dist)

{ // aka phi

   return dist.shape();

}

template <class RealType, class Policy>

inline RealType standard_deviation(const inverse_gaussian_distribution<RealType, Policy>& dist)

{

  BOOST_MATH_STD_USING

  RealType scale = dist.scale();

  RealType mean = dist.mean();

  RealType result = sqrt(mean * mean * mean / scale);

  return result;

}

template <class RealType, class Policy>

inline RealType mode(const inverse_gaussian_distribution<RealType, Policy>& dist)

{

  BOOST_MATH_STD_USING

  RealType scale = dist.scale();

  RealType  mean = dist.mean();

  RealType result = mean * (sqrt(1 + (9 * mean * mean)/(4 * scale * scale)) 

      - 3 * mean / (2 * scale));

  return result;

}

template <class RealType, class Policy>

inline RealType skewness(const inverse_gaussian_distribution<RealType, Policy>& dist)

{

  BOOST_MATH_STD_USING

  RealType scale = dist.scale();

  RealType  mean = dist.mean();

  RealType result = 3 * sqrt(mean/scale);

  return result;

}

template <class RealType, class Policy>

inline RealType kurtosis(const inverse_gaussian_distribution<RealType, Policy>& dist)

{

  RealType scale = dist.scale();

  RealType  mean = dist.mean();

  RealType result = 15 * mean / scale -3;

  return result;

}

template <class RealType, class Policy>

inline RealType kurtosis_excess(const inverse_gaussian_distribution<RealType, Policy>& dist)

{

  RealType scale = dist.scale();

  RealType  mean = dist.mean();

  RealType result = 15 * mean / scale;

  return result;

}

} // namespace math

} // namespace boost

// This include must be at the end, *after* the accessors

// for this distribution have been defined, in order to

// keep compilers that support two-phase lookup happy.

#include <boost/math/distributions/detail/derived_accessors.hpp>

#endif // BOOST_STATS_INVERSE_GAUSSIAN_HPP