201 lines · plain
1// Copyright John Maddock 2007.2// Copyright Paul A. Bristow 2007.3 4// Use, modification and distribution are subject to the5// Boost Software License, Version 1.0. (See accompanying file6// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_STATS_FIND_SCALE_HPP9#define BOOST_STATS_FIND_SCALE_HPP10 11#include <boost/math/distributions/fwd.hpp> // for all distribution signatures.12#include <boost/math/distributions/complement.hpp>13#include <boost/math/policies/policy.hpp>14// using boost::math::policies::policy;15#include <boost/math/tools/traits.hpp>16#include <boost/math/tools/assert.hpp>17#include <boost/math/special_functions/fpclassify.hpp>18#include <boost/math/policies/error_handling.hpp>19// using boost::math::complement; // will be needed by users who want complement,20// but NOT placed here to avoid putting it in global scope.21 22namespace boost23{24 namespace math25 {26 // Function to find location of random variable z27 // to give probability p (given scale)28 // Applies to normal, lognormal, extreme value, Cauchy, (and symmetrical triangular),29 // distributions that have scale.30 // BOOST_STATIC_ASSERTs, see below, are used to enforce this.31 32 template <class Dist, class Policy>33 inline34 typename Dist::value_type find_scale( // For example, normal mean.35 typename Dist::value_type z, // location of random variable z to give probability, P(X > z) == p.36 // For example, a nominal minimum acceptable weight z, so that p * 100 % are > z37 typename Dist::value_type p, // probability value desired at x, say 0.95 for 95% > z.38 typename Dist::value_type location, // location parameter, for example, normal distribution mean.39 const Policy& pol 40 )41 {42 static_assert(::boost::math::tools::is_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a distribution."); 43 static_assert(::boost::math::tools::is_scaled_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a scaled distribution."); 44 static const char* function = "boost::math::find_scale<Dist, Policy>(%1%, %1%, %1%, Policy)";45 46 if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1))47 {48 return policies::raise_domain_error<typename Dist::value_type>(49 function, "Probability parameter was %1%, but must be >= 0 and <= 1!", p, pol);50 }51 if(!(boost::math::isfinite)(z))52 {53 return policies::raise_domain_error<typename Dist::value_type>(54 function, "find_scale z parameter was %1%, but must be finite!", z, pol);55 }56 if(!(boost::math::isfinite)(location))57 {58 return policies::raise_domain_error<typename Dist::value_type>(59 function, "find_scale location parameter was %1%, but must be finite!", location, pol);60 }61 62 //cout << "z " << z << ", p " << p << ", quantile(Dist(), p) "63 //<< quantile(Dist(), p) << ", z - mean " << z - location 64 //<<", sd " << (z - location) / quantile(Dist(), p) << endl;65 66 //quantile(N01, 0.001) -3.0902367 //quantile(N01, 0.01) -2.3263568 //quantile(N01, 0.05) -1.6448569 //quantile(N01, 0.333333) -0.43072870 //quantile(N01, 0.5) 0 71 //quantile(N01, 0.666667) 0.43072872 //quantile(N01, 0.9) 1.2815573 //quantile(N01, 0.95) 1.6448574 //quantile(N01, 0.99) 2.3263575 //quantile(N01, 0.999) 3.0902376 77 typename Dist::value_type result = 78 (z - location) // difference between desired x and current location.79 / quantile(Dist(), p); // standard distribution.80 81 if (result <= 0)82 { // If policy isn't to throw, return the scale <= 0.83 policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE84 }85 return result;86 } // template <class Dist, class Policy> find_scale87 88 template <class Dist>89 inline // with default policy.90 typename Dist::value_type find_scale( // For example, normal mean.91 typename Dist::value_type z, // location of random variable z to give probability, P(X > z) == p.92 // For example, a nominal minimum acceptable z, so that p * 100 % are > z93 typename Dist::value_type p, // probability value desired at x, say 0.95 for 95% > z.94 typename Dist::value_type location) // location parameter, for example, mean.95 { // Forward to find_scale using the default policy.96 return (find_scale<Dist>(z, p, location, policies::policy<>()));97 } // find_scale98 99 template <class Dist, class Real1, class Real2, class Real3, class Policy>100 inline typename Dist::value_type find_scale(101 complemented4_type<Real1, Real2, Real3, Policy> const& c)102 {103 //cout << "cparam1 q " << c.param1 // q104 // << ", c.dist z " << c.dist // z105 // << ", c.param2 l " << c.param2 // l106 // << ", quantile (Dist(), c.param1 = q) "107 // << quantile(Dist(), c.param1) //q108 // << endl;109 110 static_assert(::boost::math::tools::is_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a distribution."); 111 static_assert(::boost::math::tools::is_scaled_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a scaled distribution."); 112 static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))";113 114 // Checks on arguments, as not complemented version,115 // Explicit policy.116 typename Dist::value_type q = c.param1;117 if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1))118 {119 return policies::raise_domain_error<typename Dist::value_type>(120 function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, c.param3);121 }122 typename Dist::value_type z = c.dist;123 if(!(boost::math::isfinite)(z))124 {125 return policies::raise_domain_error<typename Dist::value_type>(126 function, "find_scale z parameter was %1%, but must be finite!", z, c.param3);127 }128 typename Dist::value_type location = c.param2;129 if(!(boost::math::isfinite)(location))130 {131 return policies::raise_domain_error<typename Dist::value_type>(132 function, "find_scale location parameter was %1%, but must be finite!", location, c.param3);133 }134 135 typename Dist::value_type result = 136 (c.dist - c.param2) // difference between desired x and current location.137 / quantile(complement(Dist(), c.param1));138 // ( z - location) / (quantile(complement(Dist(), q)) 139 if (result <= 0)140 { // If policy isn't to throw, return the scale <= 0.141 policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE142 }143 return result;144 } // template <class Dist, class Policy, class Real1, class Real2, class Real3> typename Dist::value_type find_scale145 146 // So the user can start from the complement q = (1 - p) of the probability p,147 // for example, s = find_scale<normal>(complement(z, q, l));148 149 template <class Dist, class Real1, class Real2, class Real3>150 inline typename Dist::value_type find_scale(151 complemented3_type<Real1, Real2, Real3> const& c)152 {153 //cout << "cparam1 q " << c.param1 // q154 // << ", c.dist z " << c.dist // z155 // << ", c.param2 l " << c.param2 // l156 // << ", quantile (Dist(), c.param1 = q) "157 // << quantile(Dist(), c.param1) //q158 // << endl;159 160 static_assert(::boost::math::tools::is_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a distribution."); 161 static_assert(::boost::math::tools::is_scaled_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a scaled distribution."); 162 static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))";163 164 // Checks on arguments, as not complemented version,165 // default policy policies::policy<>().166 typename Dist::value_type q = c.param1;167 if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1))168 {169 return policies::raise_domain_error<typename Dist::value_type>(170 function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, policies::policy<>());171 }172 typename Dist::value_type z = c.dist;173 if(!(boost::math::isfinite)(z))174 {175 return policies::raise_domain_error<typename Dist::value_type>(176 function, "find_scale z parameter was %1%, but must be finite!", z, policies::policy<>());177 }178 typename Dist::value_type location = c.param2;179 if(!(boost::math::isfinite)(location))180 {181 return policies::raise_domain_error<typename Dist::value_type>(182 function, "find_scale location parameter was %1%, but must be finite!", location, policies::policy<>());183 }184 185 typename Dist::value_type result = 186 (z - location) // difference between desired x and current location.187 / quantile(complement(Dist(), q));188 // ( z - location) / (quantile(complement(Dist(), q)) 189 if (result <= 0)190 { // If policy isn't to throw, return the scale <= 0.191 policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", // LCOV_EXCL_LINE192 result, policies::policy<>()); // This is only the default policy - also Want a version with Policy here. LCOV_EXCL_LINE193 }194 return result;195 } // template <class Dist, class Real1, class Real2, class Real3> typename Dist::value_type find_scale196 197 } // namespace boost198} // namespace math199 200#endif // BOOST_STATS_FIND_SCALE_HPP201