Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r76935 - in sandbox/math/boost/math: distributions special_functions
From: ben_at_[hidden]
Date: 2012-02-07 12:12:16


Author: bensobotta
Date: 2012-02-07 12:12:14 EST (Tue, 07 Feb 2012)
New Revision: 76935
URL: http://svn.boost.org/trac/boost/changeset/76935

Log:
Many constants have been substituted with the newly added in trunk.
Did some changes to the mode computation routine, mostly the starting
guess.

Text files modified:
   sandbox/math/boost/math/distributions/skew_normal.hpp | 252 ++++++++++++++++++++++++++-------------
   sandbox/math/boost/math/special_functions/owens_t.hpp | 30 +++-
   2 files changed, 186 insertions(+), 96 deletions(-)

Modified: sandbox/math/boost/math/distributions/skew_normal.hpp
==============================================================================
--- sandbox/math/boost/math/distributions/skew_normal.hpp (original)
+++ sandbox/math/boost/math/distributions/skew_normal.hpp 2012-02-07 12:12:14 EST (Tue, 07 Feb 2012)
@@ -21,7 +21,8 @@
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/tools/tuple.hpp>
 #include <boost/math/tools/roots.hpp> // Newton-Raphson
-//#include <boost/math/distributions/detail/generic_mode.hpp> // pdf max finder.
+#include <boost/assert.hpp>
+#include <boost/math/distributions/detail/generic_mode.hpp> // pdf max finder.
 
 #include <utility>
 #include <algorithm> // std::lower_bound, std::distance
@@ -200,7 +201,7 @@
 
     normal_distribution<RealType, Policy> std_normal;
 
- result = cdf(std_normal, transformed_x) - owens_t(transformed_x, shape)*2;
+ result = cdf(std_normal, transformed_x) - owens_t(transformed_x, shape)*static_cast<RealType>(2);
 
     return result;
   } // cdf
@@ -243,7 +244,7 @@
 
     normal_distribution<RealType, Policy> std_normal;
 
- result = cdf(complement(std_normal, transformed_x)) + owens_t(transformed_x, shape)*2;
+ result = cdf(complement(std_normal, transformed_x)) + owens_t(transformed_x, shape)*static_cast<RealType>(2);
     return result;
   } // cdf complement
 
@@ -271,12 +272,10 @@
     BOOST_MATH_STD_USING // for ADL of std functions
 
     using namespace boost::math::constants;
- // TODO replace by root_two_div_pi<RealType>() when in boost-trunk
 
     const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape());
 
- return dist.location() + dist.scale()*delta*root_two<RealType>()/root_pi<RealType>();
- //return dist.location() + dist.scale() * delta * root_two_div_pi<RealType>();
+ return dist.location() + dist.scale() * delta * root_two_div_pi<RealType>();
   }
 
   template <class RealType, class Policy>
@@ -285,21 +284,135 @@
     BOOST_MATH_STD_USING // for ADL of std functions
 
     using namespace boost::math::constants;
- // TODO replace by two_div_pi<RealType>() when in boost-trunk
 
     const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape());
- RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-static_cast<RealType>(2)/pi<RealType>()*delta*delta);
- //RealType variance = dist.scale()*dist.scale()*(1-two_div_pi<RealType>() *delta*delta);
+ RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()*delta*delta);
     return variance;
   }
 
- /*
- * TODO No closed expression, so use f'(x) = 0.
- */
-
   namespace detail
   {
-
+ /*
+ TODO No closed expression for mode, so use max of pdf.
+ */
+
+ template <class RealType, class Policy>
+ inline RealType mode_fallback(const skew_normal_distribution<RealType, Policy>& dist)
+ { // mode.
+ static const char* function = "mode(skew_normal_distribution<%1%> const&)";
+ const RealType scale = dist.scale();
+ const RealType location = dist.location();
+ const RealType shape = dist.shape();
+
+ RealType result;
+ if(!detail::check_scale(
+ function,
+ scale, &result, Policy())
+ ||
+ !detail::check_skew_normal_shape(
+ function,
+ shape,
+ &result,
+ Policy()))
+ return result;
+
+ if( shape == 0 )
+ {
+ return location;
+ }
+
+ if( shape < 0 )
+ {
+ skew_normal_distribution<RealType, Policy> D(0, 1, -shape);
+ result = mode_fallback(D);
+ result = location-scale*result;
+ return result;
+ }
+
+ BOOST_MATH_STD_USING
+
+ // 21 elements
+ static const RealType shapes[] = {
+ 0.0,
+ 1.000000000000000e-004,
+ 2.069138081114790e-004,
+ 4.281332398719396e-004,
+ 8.858667904100824e-004,
+ 1.832980710832436e-003,
+ 3.792690190732250e-003,
+ 7.847599703514606e-003,
+ 1.623776739188722e-002,
+ 3.359818286283781e-002,
+ 6.951927961775606e-002,
+ 1.438449888287663e-001,
+ 2.976351441631319e-001,
+ 6.158482110660261e-001,
+ 1.274274985703135e+000,
+ 2.636650898730361e+000,
+ 5.455594781168514e+000,
+ 1.128837891684688e+001,
+ 2.335721469090121e+001,
+ 4.832930238571753e+001,
+ 1.000000000000000e+002};
+
+ // 21 elements
+ static const RealType guess[] = {
+ 0.0,
+ 5.000050000525391e-005,
+ 1.500015000148736e-004,
+ 3.500035000350010e-004,
+ 7.500075000752560e-004,
+ 1.450014500145258e-003,
+ 3.050030500305390e-003,
+ 6.250062500624765e-003,
+ 1.295012950129504e-002,
+ 2.675026750267495e-002,
+ 5.525055250552491e-002,
+ 1.132511325113255e-001,
+ 2.249522495224952e-001,
+ 3.992539925399257e-001,
+ 5.353553535535358e-001,
+ 4.954549545495457e-001,
+ 3.524535245352451e-001,
+ 2.182521825218249e-001,
+ 1.256512565125654e-001,
+ 6.945069450694508e-002,
+ 3.735037350373460e-002
+ };
+
+ const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape);
+
+ typedef typename std::iterator_traits<RealType*>::difference_type diff_type;
+
+ const diff_type d = std::distance(shapes, result_ptr);
+
+ BOOST_ASSERT(d > static_cast<diff_type>(0));
+
+ // refine
+ if(d < static_cast<diff_type>(21)) // shape smaller 100
+ {
+ result = guess[d-static_cast<diff_type>(1)]
+ + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)])
+ * (shape-shapes[d-static_cast<diff_type>(1)]);
+ }
+ else // shape greater 100
+ {
+ result = 1e-4;
+ }
+
+ skew_normal_distribution<RealType, Policy> helper(0, 1, shape);
+
+ result = detail::generic_find_mode_01(helper, result, function);
+
+ result = result*scale + location;
+
+ return result;
+ } // mode_fallback
+
+
+ /*
+ * TODO No closed expression for mode, so use f'(x) = 0
+ */
     template <class RealType, class Policy>
     struct skew_normal_mode_functor
     {
@@ -323,7 +436,7 @@
     private:
       const boost::math::skew_normal_distribution<RealType, Policy> distribution;
     };
-
+
   } // namespace detail
   
   template <class RealType, class Policy>
@@ -348,11 +461,6 @@
       return location;
     }
 
- if((boost::math::isinf)(shape))
- {
- return location;
- }
-
     if( shape < 0 )
     {
       skew_normal_distribution<RealType, Policy> D(0, 1, -shape);
@@ -361,8 +469,9 @@
       return result;
     }
 
- // 20 elements
+ // 21 elements
     static const RealType shapes[] = {
+ 0.0,
       1.000000000000000e-004,
       2.069138081114790e-004,
       4.281332398719396e-004,
@@ -386,6 +495,7 @@
 
     // 21 elements
     static const RealType guess[] = {
+ 0.0,
       5.000050000525391e-005,
       1.500015000148736e-004,
       3.500035000350010e-004,
@@ -405,18 +515,35 @@
       2.182521825218249e-001,
       1.256512565125654e-001,
       6.945069450694508e-002,
- 3.735037350373460e-002,
- 1.0e-004
+ 3.735037350373460e-002
     };
 
- const RealType* result_ptr = std::lower_bound(shapes, shapes+20, shape);
+ const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape);
 
- result = guess[std::distance(shapes, result_ptr)];
+ typedef typename std::iterator_traits<RealType*>::difference_type diff_type;
+
+ const diff_type d = std::distance(shapes, result_ptr);
+
+ BOOST_ASSERT(d > static_cast<diff_type>(0));
 
     // TODO: make the search bounds smarter, depending on the shape parameter
- const RealType search_min = 0; // below zero was caught above
- const RealType search_max = 0.55; // will never go above 0.55
+ RealType search_min = 0; // below zero was caught above
+ RealType search_max = 0.55; // will never go above 0.55
 
+ // refine
+ if(d < static_cast<diff_type>(21)) // shape smaller 100
+ {
+ // it is safe to assume that d > 0, because shape==0.0 is caught earlier
+ result = guess[d-static_cast<diff_type>(1)]
+ + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)])
+ * (shape-shapes[d-static_cast<diff_type>(1)]);
+ }
+ else // shape greater 100
+ {
+ result = 1e-4;
+ search_max = guess[19]; // set 19 instead of 20 to have a safety margin because the table may not be exact @ shape=100
+ }
+
     const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
     boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
 
@@ -424,93 +551,46 @@
 
     result = tools::newton_raphson_iterate(detail::skew_normal_mode_functor<RealType, Policy>(helper), result,
       search_min, search_max, get_digits, m);
-
+
     result = result*scale + location;
 
     return result;
-}
-
- /*
- TODO No closed expression, so use max of pdf.
- */
-
- /*
-
- template <class RealType, class Policy>
- inline RealType mode(const skew_normal_distribution<RealType, Policy>& dist)
- { // mode.
- static const char* function = "mode(skew_normal_distribution<%1%> const&)";
- RealType scale = dist.scale();
- RealType shape = dist.shape();
- RealType r;
- if(!detail::check_scale(
- function,
- scale, &r, Policy())
- ||
- !detail::check_finite(
- function,
- shape,
- &r,
- Policy()))
- return (RealType)r;
-
- BOOST_MATH_STD_USING
-
- //RealType m = v < 3 ? 0 : detail::mean(v, l, Policy());
- //RealType var = v < 4 ? 1 : detail::variance(v, l, Policy());
-
- RealType guess = (scale < 3) ? 0 : mean(dist);
- RealType step = (scale < 4) ? 0 : variance(dist); //
-
- cout << "guess= " << guess << ", step = " << step << endl;
-
- // TODO problem here in finding lower and upper bounds.
- // Guess looks reasonable but What should step be?
-
- return detail::generic_find_mode(
- dist,
- guess,
- function);
- } // mode
+ }
   
- */
 
+
   template <class RealType, class Policy>
   inline RealType skewness(const skew_normal_distribution<RealType, Policy>& dist)
   {
     BOOST_MATH_STD_USING // for ADL of std functions
     using namespace boost::math::constants;
- // TODO replace by two_div_pi<RealType>() when in boost-trunk
- // static const RealType two_over_pi = two_div_pi<RealType>(); // = 2 / pi
- static const RealType two_over_pi = static_cast<RealType>(2)/pi<RealType>();
- static const RealType sqrt_two_over_pi = root_two<RealType>()/root_pi<RealType>();
+
     static const RealType factor = four_minus_pi<RealType>()/static_cast<RealType>(2);
     const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape());
 
- return factor * pow(sqrt_two_over_pi * delta, 3) /
- pow(static_cast<RealType>(1)-two_over_pi*delta*delta, static_cast<RealType>(1.5));
+ return factor * pow(root_two_div_pi<RealType>() * delta, 3) /
+ pow(static_cast<RealType>(1)-two_div_pi<RealType>()*delta*delta, static_cast<RealType>(1.5));
   }
 
   template <class RealType, class Policy>
   inline RealType kurtosis(const skew_normal_distribution<RealType, Policy>& dist)
   {
- return kurtosis_excess(dist)+3;
+ return kurtosis_excess(dist)+static_cast<RealType>(3);
   }
 
   template <class RealType, class Policy>
   inline RealType kurtosis_excess(const skew_normal_distribution<RealType, Policy>& dist)
   {
     BOOST_MATH_STD_USING // for ADL of std functions
- // TODO replace by two_div_pi<RealType>() etc when in boost-trunk
+
     using namespace boost::math::constants;
 
- static const RealType two_over_pi = static_cast<RealType>(2)/pi<RealType>();
     static const RealType factor = pi_minus_three<RealType>()*static_cast<RealType>(2);
 
     const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape());
 
- const RealType x = static_cast<RealType>(1)-two_over_pi*delta*delta;
- const RealType y = two_over_pi * delta * delta;
+ const RealType x = static_cast<RealType>(1)-two_div_pi<RealType>()*delta*delta;
+ const RealType y = two_div_pi<RealType>() * delta * delta;
 
     return factor * y*y / (x*x);
   }
@@ -569,7 +649,9 @@
       const RealType skew = skewness(dist);
       const RealType exk = kurtosis_excess(dist);
 
- x = x + (x*x-1)*skew/6 + x*(x*x-3)*exk/24 - x*(2*x*x-5)*skew*skew/36;
+ x = x + (x*x-static_cast<RealType>(1))*skew/static_cast<RealType>(6)
+ + x*(x*x-static_cast<RealType>(3))*exk/static_cast<RealType>(24)
+ - x*(static_cast<RealType>(2)*x*x-static_cast<RealType>(5))*skew*skew/static_cast<RealType>(36);
     } // if(shape != 0)
 
     result = standard_deviation(dist)*x+mean(dist);

Modified: sandbox/math/boost/math/special_functions/owens_t.hpp
==============================================================================
--- sandbox/math/boost/math/special_functions/owens_t.hpp (original)
+++ sandbox/math/boost/math/special_functions/owens_t.hpp 2012-02-07 12:12:14 EST (Tue, 07 Feb 2012)
@@ -36,7 +36,7 @@
          inline RealType owens_t_znorm1(const RealType x)
          {
             using namespace boost::math::constants;
- return erf(x/root_two<RealType>())*half<RealType>();
+ return erf(x*one_div_root_two<RealType>())*half<RealType>();
          } // RealType owens_t_znorm1(const RealType x)
 
          // owens_t_znorm2(x) = P(x<=Z<oo) with Z being normally distributed.
@@ -44,7 +44,7 @@
          inline RealType owens_t_znorm2(const RealType x)
          {
             using namespace boost::math::constants;
- return erfc(x/root_two<RealType>())*half<RealType>();
+ return erfc(x*one_div_root_two<RealType>())*half<RealType>();
          } // RealType owens_t_znorm2(const RealType x)
 
          // Auxiliary function, it computes an array key that is used to determine
@@ -118,11 +118,11 @@
 
             unsigned short j=1;
             RealType jj = 1;
- RealType aj = a / (static_cast<RealType>(2)*pi<RealType>());
+ RealType aj = a * one_div_two_pi<RealType>();
             RealType dj = expm1( hs );
             RealType gj = hs*dhs;
 
- RealType val = atan( a ) / (static_cast<RealType>(2)*pi<RealType>());
+ RealType val = atan( a ) * one_div_two_pi<RealType>();
 
             while( true )
             {
@@ -155,7 +155,7 @@
 
             unsigned short ii = 1;
             RealType val = 0;
- RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
+ RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
             RealType z = owens_t_znorm1(ah)/h;
 
             while( true )
@@ -163,7 +163,7 @@
                val += z;
                if( maxii <= ii )
                {
- val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();
+ val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
                   break;
                } // if( maxii <= ii )
                z = y * ( vi - static_cast<RealType>(ii) * z );
@@ -202,7 +202,7 @@
 
             RealType ii = 1;
             unsigned short i = 0;
- RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
+ RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
             RealType zi = owens_t_znorm1(ah)/h;
             RealType val = 0;
 
@@ -212,7 +212,7 @@
                val += zi*c2[i];
                if( m <= i ) // if( m < i+1 )
                {
- val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();
+ val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
                   break;
                } // if( m < i )
                zi = y * (ii*zi - vi);
@@ -236,7 +236,7 @@
             const RealType as = -a*a;
 
             unsigned short ii = 1;
- RealType ai = a * exp( -hs*(static_cast<RealType>(1)-as)*half<RealType>() ) / (static_cast<RealType>(2)*pi<RealType>());
+ RealType ai = a * exp( -hs*(static_cast<RealType>(1)-as)*half<RealType>() ) * one_div_two_pi<RealType>();
             RealType yi = 1;
             RealType val = 0;
 
@@ -258,6 +258,14 @@
          inline RealType owens_t_T5(const RealType h, const RealType a, const unsigned short m)
          {
             BOOST_MATH_STD_USING
+ /*
+ NOTICE:
+ - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
+ polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
+ quadrature, because T5(h,a,m) contains only x^2 terms.
+ - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
+ of 1/(2*pi) according to T5(h,a,m).
+ */
             static const RealType pts[] = {0.35082039676451715489E-02,
                0.31279042338030753740E-01, 0.85266826283219451090E-01,
                0.16245071730812277011, 0.25851196049125434828,
@@ -301,7 +309,7 @@
             RealType val = normh * ( static_cast<RealType>(1) - normh ) * half<RealType>();
 
             if( r != 0 )
- val -= r * exp( -y*h*h*half<RealType>()/r ) / (static_cast<RealType>(2)*pi<RealType>());
+ val -= r * exp( -y*h*h*half<RealType>()/r ) * one_div_two_pi<RealType>();
 
             return val;
          } // RealType owens_t_T6(const RealType h, const RealType a, const unsigned short m)
@@ -402,7 +410,7 @@
       {
          typedef typename tools::promote_args<T1, T2>::type result_type;
          typedef typename policies::evaluation<result_type, Policy>::type value_type;
- return policies::checked_narrowing_cast<result_type, Policy>(detail::owens_t(static_cast<value_type>(h), static_cast<value_type>(a), pol), "boost::math::ellint_2<%1%>(%1%,%1%)");
+ return policies::checked_narrowing_cast<result_type, Policy>(detail::owens_t(static_cast<value_type>(h), static_cast<value_type>(a), pol), "boost::math::owens_t<%1%>(%1%,%1%)");
       }
 
       template <class T1, class T2>


Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk