Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r66766 - trunk/boost/math/distributions
From: pbristow_at_[hidden]
Date: 2010-11-26 05:29:01


Author: pbristow
Date: 2010-11-26 05:28:58 EST (Fri, 26 Nov 2010)
New Revision: 66766
URL: http://svn.boost.org/trac/boost/changeset/66766

Log:
removed diagnostics
Text files modified:
   trunk/boost/math/distributions/inverse_gaussian.hpp | 80 +---------------------------------------
   1 files changed, 2 insertions(+), 78 deletions(-)

Modified: trunk/boost/math/distributions/inverse_gaussian.hpp
==============================================================================
--- trunk/boost/math/distributions/inverse_gaussian.hpp (original)
+++ trunk/boost/math/distributions/inverse_gaussian.hpp 2010-11-26 05:28:58 EST (Fri, 26 Nov 2010)
@@ -198,23 +198,8 @@
    RealType expfactor = exp(2 * scale / mean);
    RealType n3 = - sqrt(scale / x);
    n3 *= (x / mean) + 1;
- // cout << "((x / mean) +1) = " << n3 << endl;
    RealType n4 = cdf(n01, n3);
- //cout << "phi1 = " << n1 << ", exp(2 * scale / mean) = " << n2 << ", exp * phi2 = " << n4 * n2 << endl;
-
    result = n1 + expfactor * n4;
-
- if(false)
- { // Output some diagnostic values.
- cout <<"_\n cdf===========================" << endl;
- cout << "sqrt(scale / x)*((x / mean) -1) = " << n0 << endl;
- cout << "cdf(n01, n1) = " << n1 << endl;
- cout << "exp(2 * scale / mean) = " << expfactor << endl;
- cout << " - sqrt(scale / x)*((x / mean) +1) = " << n3 << endl;
- cout << "cdf(n01, - sqrt(scale / x)*((x / mean) +1)) = " << n4 << endl;
- cout << "exp * cdf_2 = " << n4 * expfactor << endl;
- cout << "cdf_1 + exp * cdf_2 = " << result << endl;
- }
    return result;
 } // cdf
 
@@ -231,10 +216,6 @@
     RealType c = cdf(distribution, x);
     RealType fx = c - prob; // Difference cdf - value - to minimize.
     RealType dx = pdf(distribution, x); // pdf is 1st derivative.
- if(false)
- {
- cout << "cdf(dist, " << x << ") = " << c << ", diff = " << fx << ", dx " << dx << endl;
- }
     // return both function evaluation difference f(x) and 1st derivative f'(x).
     return std::tr1::make_tuple(fx, dx);
   }
@@ -255,18 +236,6 @@
     RealType c = cdf(complement(distribution, x));
     RealType fx = c - prob; // Difference cdf - value - to minimize.
     RealType dx = -pdf(distribution, x); // pdf is 1st derivative.
- if(true)
- {
- std::streamsize precision = cout.precision(); // Save
- if (false)
- {
- cout << setprecision(numeric_limits<RealType>::max_digits10)
- << "cdf((complement(dist, " << x << ")) = " << c
- << setprecision(4)
- << ", diff = " << fx << ", dx " << dx << endl;
- cout.precision(precision); // restore
- }
- }
     // return both function evaluation difference f(x) and 1st derivative f'(x).
     return std::tr1::make_tuple(fx, dx);
   }
@@ -299,20 +268,10 @@
       // 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.
- //x = qnorm(p, 0, 1, true, false);
- //x /= sqrt(phi);
- //x = x - 1. / (2 * phi);
- //x = mu * exp(x);
- // x = mu * exp(qnorm(p, 0, 1, true, false) / sqrt(phi) - 1/(2 * phi));
+
       normal_distribution<RealType, no_overthrow_policy> n01;
       x = mu * exp(quantile(n01, p) / sqrt(phi) - 1/(2 * phi));
- // Might add about 0.006 to this to get closer?
- //RealType pi = 3.1459;
- //cout << "1 / sqrt(8 * pi * phi) = " << 1 / sqrt(8 * pi * phi) << endl;
- // Bagshaw guess is:
- // RealType U = quantile(n01, p); // U <- qnorm (p)
- // RealType r1 = 1 + U / sqrt (phi) + U * U / (2 * phi) + U * U * U / (8 * phi * sqrt(phi));
- }
+ }
     else
     { // phi < 2 so much less symmetrical with long tail,
       // so use gamma distribution as an approximation.
@@ -327,15 +286,12 @@
       // 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);
- //cout << "quantile(complement(g, p)) = " << qg << ", qgamma(1.- p, 0.5, 1.0, true, false); = " << qg1 << endl; // 49.014664823030209 1.#INF
       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);
- // cout << "quantile((g, p)) = " << q << endl;// 49.014664823030209
- //<< ", qgamma(1.- p, 0.5, 1.0); = " << x << endl; // 1.#INF
        // 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));
@@ -382,19 +338,9 @@
   // 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.
- if(false)
- {
- cout << "Probability " << p << ", guess " << guess
- << ", min " << min << ", max " << max
- //<< ", std::numeric_limits<" << typeid(RealType).name() << ">::digits = " << digits
- << ", accuracy " << get_digits << " bits."
- << ", max iterations set by policy " << m
- << endl;
- }
   using boost::math::tools::newton_raphson_iterate;
   result =
     newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType>(dist, p), guess, min, max, get_digits, m);
- //cout << m << " iterations." << endl;
    return result;
 } // quantile
 
@@ -443,18 +389,6 @@
    RealType n6 = cdf(complement(n01, +sqrt(scale/x) * ((x /mean) + 1)));
    RealType n4 = cdf(n01, n3); // =
    result = cdf_1 - expfactor * n6;
- if(false)
- {
- cout <<"_\n cdf(complement ===========================" << endl;
- cout << "sqrt(scale / x)*((x / mean) -1) = " << n0 << endl;
- cout << "cdf(complement(n01, n1)) = " << cdf_1 << endl;
- cout << "-sqrt(scale / x) * ((x / mean) +1) = " << n3 << endl;
- cout << "exp(2 * scale / mean) = " << expfactor << endl;
- cout << "cdf(complement(n01, +sqrt(scale/x) * ((x /mean) + 1))) = " << n6 << endl;
- cout << "cdf((n01, ) exp(2 * scale / mean) * (x / mean) + 1) = " << n4 << endl;
- cout << "exp * cdf_2 = " << result << endl;
- }
- //cout << "cdf(complement) result = " << result << endl;
    return result;
 } // cdf complement
 
@@ -485,19 +419,9 @@
   // 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>();
- if(false)
- {
- cout << "Probability " << q << ", guess at x = " << guess
- //<< ", min " << min << ", max " << max
- ////<< ", std::numeric_limits<" << typeid(RealType).name() << ">::digits = " << digits
- // << ", accuracy " << get_digits << " bits."
- // << ", max iterations set by policy " << m
- << endl;
- }
   using boost::math::tools::newton_raphson_iterate;
   result =
     newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType>(c.dist, q), guess, min, max, get_digits, m);
- //cout << m << " iterations." << endl;
    return result;
 } // quantile
 


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