math/distributions/find_scale.hpp

0001 //  Copyright John Maddock 2007.
0002 //  Copyright Paul A. Bristow 2007.
0003
0004 //  Use, modification and distribution are subject to the
0005 //  Boost Software License, Version 1.0. (See accompanying file
0006 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_STATS_FIND_SCALE_HPP
0009 #define BOOST_STATS_FIND_SCALE_HPP
0010
0011 #include <boost/math/distributions/fwd.hpp> // for all distribution signatures.
0012 #include <boost/math/distributions/complement.hpp>
0013 #include <boost/math/policies/policy.hpp>
0014 // using boost::math::policies::policy;
0015 #include <boost/math/tools/traits.hpp>
0016 #include <boost/math/tools/assert.hpp>
0017 #include <boost/math/special_functions/fpclassify.hpp>
0018 #include <boost/math/policies/error_handling.hpp>
0019 // using boost::math::complement; // will be needed by users who want complement,
0020 // but NOT placed here to avoid putting it in global scope.
0021
0022 namespace boost
0023 {
0024   namespace math
0025   {
0026     // Function to find location of random variable z
0027     // to give probability p (given scale)
0028     // Applies to normal, lognormal, extreme value, Cauchy, (and symmetrical triangular),
0029     // distributions that have scale.
0030     // BOOST_STATIC_ASSERTs, see below, are used to enforce this.
0031
0032     template <class Dist, class Policy>
0033     inline
0034       typename Dist::value_type find_scale( // For example, normal mean.
0035       typename Dist::value_type z, // location of random variable z to give probability, P(X > z) == p.
0036       // For example, a nominal minimum acceptable weight z, so that p * 100 % are > z
0037       typename Dist::value_type p, // probability value desired at x, say 0.95 for 95% > z.
0038       typename Dist::value_type location, // location parameter, for example, normal distribution mean.
0039       const Policy& pol
0040       )
0041     {
0042       static_assert(::boost::math::tools::is_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a distribution.");
0043       static_assert(::boost::math::tools::is_scaled_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a scaled distribution.");
0044       static const char* function = "boost::math::find_scale<Dist, Policy>(%1%, %1%, %1%, Policy)";
0045
0046       if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1))
0047       {
0048         return policies::raise_domain_error<typename Dist::value_type>(
0049           function, "Probability parameter was %1%, but must be >= 0 and <= 1!", p, pol);
0050       }
0051       if(!(boost::math::isfinite)(z))
0052       {
0053         return policies::raise_domain_error<typename Dist::value_type>(
0054           function, "find_scale z parameter was %1%, but must be finite!", z, pol);
0055       }
0056       if(!(boost::math::isfinite)(location))
0057       {
0058         return policies::raise_domain_error<typename Dist::value_type>(
0059           function, "find_scale location parameter was %1%, but must be finite!", location, pol);
0060       }
0061
0062       //cout << "z " << z << ", p " << p << ",  quantile(Dist(), p) "
0063       //<< quantile(Dist(), p) << ", z - mean " << z - location
0064       //<<", sd " << (z - location)  / quantile(Dist(), p) << endl;
0065
0066       //quantile(N01, 0.001) -3.09023
0067       //quantile(N01, 0.01) -2.32635
0068       //quantile(N01, 0.05) -1.64485
0069       //quantile(N01, 0.333333) -0.430728
0070       //quantile(N01, 0.5) 0
0071       //quantile(N01, 0.666667) 0.430728
0072       //quantile(N01, 0.9) 1.28155
0073       //quantile(N01, 0.95) 1.64485
0074       //quantile(N01, 0.99) 2.32635
0075       //quantile(N01, 0.999) 3.09023
0076
0077       typename Dist::value_type result =
0078         (z - location)  // difference between desired x and current location.
0079         / quantile(Dist(), p); // standard distribution.
0080
0081       if (result <= 0)
0082       { // If policy isn't to throw, return the scale <= 0.
0083         policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE
0084       }
0085       return result;
0086     } // template <class Dist, class Policy> find_scale
0087
0088     template <class Dist>
0089     inline // with default policy.
0090       typename Dist::value_type find_scale( // For example, normal mean.
0091       typename Dist::value_type z, // location of random variable z to give probability, P(X > z) == p.
0092       // For example, a nominal minimum acceptable z, so that p * 100 % are > z
0093       typename Dist::value_type p, // probability value desired at x, say 0.95 for 95% > z.
0094       typename Dist::value_type location) // location parameter, for example, mean.
0095     { // Forward to find_scale using the default policy.
0096       return (find_scale<Dist>(z, p, location, policies::policy<>()));
0097     } // find_scale
0098
0099     template <class Dist, class Real1, class Real2, class Real3, class Policy>
0100     inline typename Dist::value_type find_scale(
0101       complemented4_type<Real1, Real2, Real3, Policy> const& c)
0102     {
0103       //cout << "cparam1 q " << c.param1 // q
0104       //  << ", c.dist z " << c.dist // z
0105       //  << ", c.param2 l " << c.param2 // l
0106       //  << ", quantile (Dist(), c.param1 = q) "
0107       //  << quantile(Dist(), c.param1) //q
0108       //  << endl;
0109
0110       static_assert(::boost::math::tools::is_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a distribution.");
0111       static_assert(::boost::math::tools::is_scaled_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a scaled distribution.");
0112       static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))";
0113
0114       // Checks on arguments, as not complemented version,
0115       // Explicit policy.
0116       typename Dist::value_type q = c.param1;
0117       if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1))
0118       {
0119         return policies::raise_domain_error<typename Dist::value_type>(
0120           function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, c.param3);
0121       }
0122       typename Dist::value_type z = c.dist;
0123       if(!(boost::math::isfinite)(z))
0124       {
0125         return policies::raise_domain_error<typename Dist::value_type>(
0126           function, "find_scale z parameter was %1%, but must be finite!", z, c.param3);
0127       }
0128       typename Dist::value_type location = c.param2;
0129       if(!(boost::math::isfinite)(location))
0130       {
0131         return policies::raise_domain_error<typename Dist::value_type>(
0132           function, "find_scale location parameter was %1%, but must be finite!", location, c.param3);
0133       }
0134
0135       typename Dist::value_type result =
0136         (c.dist - c.param2)  // difference between desired x and current location.
0137         / quantile(complement(Dist(), c.param1));
0138       //     (  z    - location) / (quantile(complement(Dist(),  q))
0139       if (result <= 0)
0140       { // If policy isn't to throw, return the scale <= 0.
0141         policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE
0142       }
0143       return result;
0144     } // template <class Dist, class Policy, class Real1, class Real2, class Real3> typename Dist::value_type find_scale
0145
0146     // So the user can start from the complement q = (1 - p) of the probability p,
0147     // for example, s = find_scale<normal>(complement(z, q, l));
0148
0149     template <class Dist, class Real1, class Real2, class Real3>
0150     inline typename Dist::value_type find_scale(
0151       complemented3_type<Real1, Real2, Real3> const& c)
0152     {
0153       //cout << "cparam1 q " << c.param1 // q
0154       //  << ", c.dist z " << c.dist // z
0155       //  << ", c.param2 l " << c.param2 // l
0156       //  << ", quantile (Dist(), c.param1 = q) "
0157       //  << quantile(Dist(), c.param1) //q
0158       //  << endl;
0159
0160       static_assert(::boost::math::tools::is_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a distribution.");
0161       static_assert(::boost::math::tools::is_scaled_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a scaled distribution.");
0162       static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))";
0163
0164       // Checks on arguments, as not complemented version,
0165       // default policy policies::policy<>().
0166       typename Dist::value_type q = c.param1;
0167       if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1))
0168       {
0169         return policies::raise_domain_error<typename Dist::value_type>(
0170           function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, policies::policy<>());
0171       }
0172       typename Dist::value_type z = c.dist;
0173       if(!(boost::math::isfinite)(z))
0174       {
0175         return policies::raise_domain_error<typename Dist::value_type>(
0176           function, "find_scale z parameter was %1%, but must be finite!", z, policies::policy<>());
0177       }
0178       typename Dist::value_type location = c.param2;
0179       if(!(boost::math::isfinite)(location))
0180       {
0181         return policies::raise_domain_error<typename Dist::value_type>(
0182           function, "find_scale location parameter was %1%, but must be finite!", location, policies::policy<>());
0183       }
0184
0185       typename Dist::value_type result =
0186         (z - location)  // difference between desired x and current location.
0187         / quantile(complement(Dist(), q));
0188       //     (  z    - location) / (quantile(complement(Dist(),  q))
0189       if (result <= 0)
0190       { // If policy isn't to throw, return the scale <= 0.
0191         policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", // LCOV_EXCL_LINE
0192            result, policies::policy<>()); // This is only the default policy - also Want a version with Policy here.  LCOV_EXCL_LINE
0193       }
0194       return result;
0195     } // template <class Dist, class Real1, class Real2, class Real3> typename Dist::value_type find_scale
0196
0197   } // namespace boost
0198 } // namespace math
0199
0200 #endif // BOOST_STATS_FIND_SCALE_HPP