1 // Copyright John Maddock 2007. 2 // Copyright Paul A. Bristow 2007. 3 4 // Use, modification and distribution are subject to the 5 // Boost Software License, Version 1.0. (See accompanying file 6 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 7 8 #ifndef BOOST_STATS_FIND_SCALE_HPP 9 #define BOOST_STATS_FIND_SCALE_HPP 10 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 22 namespace boost 23 { 24 namespace math 25 { 26 // Function to find location of random variable z 27 // 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 inline find_scale(typename Dist::value_type z,typename Dist::value_type p,typename Dist::value_type location,const Policy & pol)34 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 > z 37 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.09023 67 //quantile(N01, 0.01) -2.32635 68 //quantile(N01, 0.05) -1.64485 69 //quantile(N01, 0.333333) -0.430728 70 //quantile(N01, 0.5) 0 71 //quantile(N01, 0.666667) 0.430728 72 //quantile(N01, 0.9) 1.28155 73 //quantile(N01, 0.95) 1.64485 74 //quantile(N01, 0.99) 2.32635 75 //quantile(N01, 0.999) 3.09023 76 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, 84 "Computed scale (%1%) is <= 0!" " Was the complement intended?", 85 result, Policy()); 86 } 87 return result; 88 } // template <class Dist, class Policy> find_scale 89 90 template <class Dist> 91 inline // with default policy. find_scale(typename Dist::value_type z,typename Dist::value_type p,typename Dist::value_type location)92 typename Dist::value_type find_scale( // For example, normal mean. 93 typename Dist::value_type z, // location of random variable z to give probability, P(X > z) == p. 94 // For example, a nominal minimum acceptable z, so that p * 100 % are > z 95 typename Dist::value_type p, // probability value desired at x, say 0.95 for 95% > z. 96 typename Dist::value_type location) // location parameter, for example, mean. 97 { // Forward to find_scale using the default policy. 98 return (find_scale<Dist>(z, p, location, policies::policy<>())); 99 } // find_scale 100 101 template <class Dist, class Real1, class Real2, class Real3, class Policy> find_scale(complemented4_type<Real1,Real2,Real3,Policy> const & c)102 inline typename Dist::value_type find_scale( 103 complemented4_type<Real1, Real2, Real3, Policy> const& c) 104 { 105 //cout << "cparam1 q " << c.param1 // q 106 // << ", c.dist z " << c.dist // z 107 // << ", c.param2 l " << c.param2 // l 108 // << ", quantile (Dist(), c.param1 = q) " 109 // << quantile(Dist(), c.param1) //q 110 // << endl; 111 112 static_assert(::boost::math::tools::is_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a distribution."); 113 static_assert(::boost::math::tools::is_scaled_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a scaled distribution."); 114 static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))"; 115 116 // Checks on arguments, as not complemented version, 117 // Explicit policy. 118 typename Dist::value_type q = c.param1; 119 if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1)) 120 { 121 return policies::raise_domain_error<typename Dist::value_type>( 122 function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, c.param3); 123 } 124 typename Dist::value_type z = c.dist; 125 if(!(boost::math::isfinite)(z)) 126 { 127 return policies::raise_domain_error<typename Dist::value_type>( 128 function, "find_scale z parameter was %1%, but must be finite!", z, c.param3); 129 } 130 typename Dist::value_type location = c.param2; 131 if(!(boost::math::isfinite)(location)) 132 { 133 return policies::raise_domain_error<typename Dist::value_type>( 134 function, "find_scale location parameter was %1%, but must be finite!", location, c.param3); 135 } 136 137 typename Dist::value_type result = 138 (c.dist - c.param2) // difference between desired x and current location. 139 / quantile(complement(Dist(), c.param1)); 140 // ( z - location) / (quantile(complement(Dist(), q)) 141 if (result <= 0) 142 { // If policy isn't to throw, return the scale <= 0. 143 policies::raise_evaluation_error<typename Dist::value_type>(function, 144 "Computed scale (%1%) is <= 0!" " Was the complement intended?", 145 result, Policy()); 146 } 147 return result; 148 } // template <class Dist, class Policy, class Real1, class Real2, class Real3> typename Dist::value_type find_scale 149 150 // So the user can start from the complement q = (1 - p) of the probability p, 151 // for example, s = find_scale<normal>(complement(z, q, l)); 152 153 template <class Dist, class Real1, class Real2, class Real3> find_scale(complemented3_type<Real1,Real2,Real3> const & c)154 inline typename Dist::value_type find_scale( 155 complemented3_type<Real1, Real2, Real3> const& c) 156 { 157 //cout << "cparam1 q " << c.param1 // q 158 // << ", c.dist z " << c.dist // z 159 // << ", c.param2 l " << c.param2 // l 160 // << ", quantile (Dist(), c.param1 = q) " 161 // << quantile(Dist(), c.param1) //q 162 // << endl; 163 164 static_assert(::boost::math::tools::is_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a distribution."); 165 static_assert(::boost::math::tools::is_scaled_distribution<Dist>::value, "The provided distribution does not meet the conceptual requirements of a scaled distribution."); 166 static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))"; 167 168 // Checks on arguments, as not complemented version, 169 // default policy policies::policy<>(). 170 typename Dist::value_type q = c.param1; 171 if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1)) 172 { 173 return policies::raise_domain_error<typename Dist::value_type>( 174 function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, policies::policy<>()); 175 } 176 typename Dist::value_type z = c.dist; 177 if(!(boost::math::isfinite)(z)) 178 { 179 return policies::raise_domain_error<typename Dist::value_type>( 180 function, "find_scale z parameter was %1%, but must be finite!", z, policies::policy<>()); 181 } 182 typename Dist::value_type location = c.param2; 183 if(!(boost::math::isfinite)(location)) 184 { 185 return policies::raise_domain_error<typename Dist::value_type>( 186 function, "find_scale location parameter was %1%, but must be finite!", location, policies::policy<>()); 187 } 188 189 typename Dist::value_type result = 190 (z - location) // difference between desired x and current location. 191 / quantile(complement(Dist(), q)); 192 // ( z - location) / (quantile(complement(Dist(), q)) 193 if (result <= 0) 194 { // If policy isn't to throw, return the scale <= 0. 195 policies::raise_evaluation_error<typename Dist::value_type>(function, 196 "Computed scale (%1%) is <= 0!" " Was the complement intended?", 197 result, policies::policy<>()); // This is only the default policy - also Want a version with Policy here. 198 } 199 return result; 200 } // template <class Dist, class Real1, class Real2, class Real3> typename Dist::value_type find_scale 201 202 } // namespace boost 203 } // namespace math 204 205 #endif // BOOST_STATS_FIND_SCALE_HPP 206