[Boost-commit] svn:boost r43222 - in sandbox/math_toolkit: boost/math/distributions libs/math/doc/sf_and_dist/graphs

From: john_at_[hidden]
Date: 2008-02-11 12:21:15

Author: johnmaddock
Date: 2008-02-11 12:21:14 EST (Mon, 11 Feb 2008)
New Revision: 43222
URL: http://svn.boost.org/trac/boost/changeset/43222

Log:
Updated non central beta cdf to accept both x and 1-x on input: allows greater accuracy when using these routines in the non central F and T distributions.
Updated non central F to take advantage of this.
Added some non central distributions to the graph generator.
Text files modified: 
   sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp    |    34 +++++++++++------                       
   sandbox/math_toolkit/boost/math/distributions/non_central_f.hpp       |    79 +++++++++++++++++++++++++++++++++++++-- 
   sandbox/math_toolkit/libs/math/doc/sf_and_dist/graphs/dist_graphs.cpp |    20 ++++++++++                              
   3 files changed, 116 insertions(+), 17 deletions(-)

Modified: sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp	(original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp	2008-02-11 12:21:14 EST (Mon, 11 Feb 2008)
@@ -30,7 +30,7 @@
       namespace detail{
 
          template <class T, class Policy>
-         T non_central_beta_p(T a, T b, T lam, T x, const Policy& pol, T init_val = 0)
+         T non_central_beta_p(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0)
          {
             BOOST_MATH_STD_USING
                using namespace boost::math;
@@ -48,9 +48,14 @@
             // Starting Poisson weight:
             T pois = gamma_p_derivative(k+1, l2, pol);
             // Starting beta term:
-            T beta = ibeta(a + k, b, x, pol);
+            T beta = x < y 
+               ? ibeta(a + k, b, x, pol)
+               : ibetac(b, a + k, y, pol);
             // recurance term:
-            T xterm = ibeta_derivative(a + k, b, x, pol) * (1 - x) / (a + b + k - 1);
+            T xterm = x < y
+               ? ibeta_derivative(a + k, b, x, pol)
+               : ibeta_derivative(b, a + k, y, pol);
+            xterm *= (1 - x) / (a + b + k - 1);
             T poisf(pois), betaf(beta), xtermf(xterm);
             T sum = init_val;
 
@@ -89,7 +94,7 @@
          }
 
          template <class T, class Policy>
-         T non_central_beta_q(T a, T b, T lam, T x, const Policy& pol, T init_val = 0)
+         T non_central_beta_q(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0)
          {
             BOOST_MATH_STD_USING
                using namespace boost::math;
@@ -107,9 +112,14 @@
             // Starting Poisson weight:
             T pois = gamma_p_derivative(k+1, l2, pol);
             // Starting beta term:
-            T beta = ibetac(a + k, b, x, pol);
+            T beta = x < y
+               ? ibetac(a + k, b, x, pol)
+               : ibeta(b, a + k, y, pol);
             // recurance term:
-            T xterm = ibeta_derivative(a + k, b, x, pol) * (1 - x) / (a + b + k - 1);
+            T xterm = x < y 
+               ? ibeta_derivative(a + k, b, x, pol)
+               : ibeta_derivative(b, a + k, y, pol);
+            xterm *= (1 - x) / (a + b + k - 1);
             T poisf(pois), betaf(beta), xtermf(xterm);
             T sum = init_val;
             //
@@ -148,7 +158,7 @@
          }
 
          template <class RealType, class Policy>
-         inline RealType non_central_beta_cdf(RealType x, RealType a, RealType b, RealType l, bool invert, const Policy&)
+         inline RealType non_central_beta_cdf(RealType x, RealType y, RealType a, RealType b, RealType l, bool invert, const Policy&)
          {
             typedef typename policies::evaluation<RealType, Policy>::type value_type;
             typedef typename policies::normalise<
@@ -162,7 +172,7 @@
 
             if(x == 0)
                return invert ? 1.0f : 0.0f;
-            if(x == 1)
+            if(y == 0)
                return invert ? 0.0f : 1.0f;
             value_type result;
             value_type c = a + b + l / 2;
@@ -177,6 +187,7 @@
                   static_cast<value_type>(b), 
                   static_cast<value_type>(l), 
                   static_cast<value_type>(x), 
+                  static_cast<value_type>(y), 
                   forwarding_policy(), 
                   static_cast<value_type>(invert ? 0 : -1));
                invert = !invert;
@@ -188,6 +199,7 @@
                   static_cast<value_type>(b), 
                   static_cast<value_type>(l), 
                   static_cast<value_type>(x), 
+                  static_cast<value_type>(y), 
                   forwarding_policy(),
                   static_cast<value_type>(invert ? -1 : 0));
             }
@@ -712,7 +724,7 @@
                Policy()))
                   return (RealType)r;
 
-         return detail::non_central_beta_cdf(x, a, b, l, false, Policy());
+         return detail::non_central_beta_cdf(x, 1 - x, a, b, l, false, Policy());
       } // cdf
 
       template <class RealType, class Policy>
@@ -746,7 +758,7 @@
                Policy()))
                   return (RealType)r;
 
-         return detail::non_central_beta_cdf(x, a, b, l, true, Policy());
+         return detail::non_central_beta_cdf(x, 1 - x, a, b, l, true, Policy());
       } // ccdf
 
       template <class RealType, class Policy>
@@ -771,5 +783,3 @@
 
 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP
 
-
-

Modified: sandbox/math_toolkit/boost/math/distributions/non_central_f.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_f.hpp	(original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_f.hpp	2008-02-11 12:21:14 EST (Mon, 11 Feb 2008)
@@ -287,24 +287,83 @@
       template <class RealType, class Policy>
       RealType cdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)
       { 
+         const char* function = "cdf(const non_central_f_distribution<%1%>&, %1%)";
+         RealType r;
+         if(!detail::check_df(
+            function,
+            dist.degrees_of_freedom1(), &r, Policy())
+               ||
+            !detail::check_df(
+               function,
+               dist.degrees_of_freedom2(), &r, Policy())
+               ||
+            !detail::check_non_centrality(
+               function,
+               dist.non_centrality(),
+               &r,
+               Policy()))
+               return r;
+         
+         if((x < 0) || !(boost::math::isfinite)(x))
+         {
+            return policies::raise_domain_error<RealType>(
+               function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());
+         }
+
          RealType alpha = dist.degrees_of_freedom1() / 2;
          RealType beta = dist.degrees_of_freedom2() / 2;
          RealType y = x * alpha / beta;
-         RealType r;
          RealType c = y / (1 + y);
-         //RealType cp = 1 / (1 + y);
-         r = cdf(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), c);
-         //RealType r2 = cdf(complement(boost::math::non_central_beta_distribution<RealType, Policy>(beta, alpha, dist.non_centrality()), cp));
+         RealType cp = 1 / (1 + y);
+         //
+         // To ensure accuracy, we pass both x and 1-x to the
+         // non-central beta cdf routine, this ensures accuracy
+         // even when we compute x to be ~ 1:
+         //
+         r = detail::non_central_beta_cdf(c, cp, alpha, beta, 
+            dist.non_centrality(), false, Policy());
          return r;
       } // cdf
 
       template <class RealType, class Policy>
       RealType cdf(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)
       { // Complemented Cumulative Distribution Function
+         const char* function = "cdf(complement(const non_central_f_distribution<%1%>&, %1%))";
+         RealType r;
+         if(!detail::check_df(
+            function,
+            c.dist.degrees_of_freedom1(), &r, Policy())
+               ||
+            !detail::check_df(
+               function,
+               c.dist.degrees_of_freedom2(), &r, Policy())
+               ||
+            !detail::check_non_centrality(
+               function,
+               c.dist.non_centrality(),
+               &r,
+               Policy()))
+               return r;
+         
+         if((c.param < 0) || !(boost::math::isfinite)(c.param))
+         {
+            return policies::raise_domain_error<RealType>(
+               function, "Random Variable parameter was %1%, but must be > 0 !", c.param, Policy());
+         }
+
          RealType alpha = c.dist.degrees_of_freedom1() / 2;
          RealType beta = c.dist.degrees_of_freedom2() / 2;
          RealType y = c.param * alpha / beta;
-         return cdf(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), y / (1 + y)));
+         RealType x = y / (1 + y);
+         RealType cx = 1 / (1 + y);
+         //
+         // To ensure accuracy, we pass both x and 1-x to the
+         // non-central beta cdf routine, this ensures accuracy
+         // even when we compute x to be ~ 1:
+         //
+         r = detail::non_central_beta_cdf(x, cx, alpha, beta, 
+            c.dist.non_centrality(), true, Policy());
+         return r;
       } // ccdf
 
       template <class RealType, class Policy>
@@ -313,6 +372,11 @@
          RealType alpha = dist.degrees_of_freedom1() / 2;
          RealType beta = dist.degrees_of_freedom2() / 2;
          RealType x = quantile(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), p);
+         if(x == 1)
+            return policies::raise_overflow_error<RealType>(
+               "quantile(const non_central_f_distribution<%1%>&, %1%)", 
+               "Result of non central F quantile is too large to represent.", 
+               Policy());
          return (x / (1 - x)) * (dist.degrees_of_freedom2() / dist.degrees_of_freedom1());
       } // quantile
 
@@ -322,6 +386,11 @@
          RealType alpha = c.dist.degrees_of_freedom1() / 2;
          RealType beta = c.dist.degrees_of_freedom2() / 2;
          RealType x = quantile(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), c.param));
+         if(x == 1)
+            return policies::raise_overflow_error<RealType>(
+               "quantile(complement(const non_central_f_distribution<%1%>&, %1%))", 
+               "Result of non central F quantile is too large to represent.", 
+               Policy());
          return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1());
       } // quantile complement.
 

Modified: sandbox/math_toolkit/libs/math/doc/sf_and_dist/graphs/dist_graphs.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/sf_and_dist/graphs/dist_graphs.cpp	(original)
+++ sandbox/math_toolkit/libs/math/doc/sf_and_dist/graphs/dist_graphs.cpp	2008-02-11 12:21:14 EST (Mon, 11 Feb 2008)
@@ -237,6 +237,26 @@
    nc_cs_plotter.add(boost::math::non_central_chi_squared(20, 100), "v=20, &#x3BB;=100");
    nc_cs_plotter.plot("Non Central Chi Squared PDF", "nccs_pdf.svg");
 
+   distribution_plotter<boost::math::non_central_beta> 
+      nc_beta_plotter;
+   nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 0), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=0");
+   nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 1), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=1");
+   nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 5), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=5");
+   nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 10), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=10");
+   nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 40), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=40");
+   nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 100), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=100");
+   nc_beta_plotter.plot("Non Central Beta PDF", "nc_beta_pdf.svg");
+
+   distribution_plotter<boost::math::non_central_f> 
+      nc_f_plotter;
+   nc_f_plotter.add(boost::math::non_central_f(10, 20, 0), "v1=10, v2=20, &#x3BB;=0");
+   nc_f_plotter.add(boost::math::non_central_f(10, 20, 1), "v1=10, v2=20, &#x3BB;=1");
+   nc_f_plotter.add(boost::math::non_central_f(10, 20, 5), "v1=10, v2=20, &#x3BB;=5");
+   nc_f_plotter.add(boost::math::non_central_f(10, 20, 10), "v1=10, v2=20, &#x3BB;=10");
+   nc_f_plotter.add(boost::math::non_central_f(10, 20, 40), "v1=10, v2=20, &#x3BB;=40");
+   nc_f_plotter.add(boost::math::non_central_f(10, 20, 100), "v1=10, v2=20, &#x3BB;=100");
+   nc_f_plotter.plot("Non Central F PDF", "nc_f_pdf.svg");
+
    distribution_plotter<boost::math::beta_distribution<> > 
       beta_plotter;
    beta_plotter.add(boost::math::beta_distribution<>(0.5, 0.5), "alpha=0.5, beta=0.5");