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/static_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 #if !defined(BOOST_NO_SFINAE) && !BOOST_WORKAROUND(__SUNPRO_CC, BOOST_TESTED_AT(0x590)) 43 BOOST_STATIC_ASSERT(::boost::math::tools::is_distribution<Dist>::value); 44 BOOST_STATIC_ASSERT(::boost::math::tools::is_scaled_distribution<Dist>::value); 45 #endif 46 static const char* function = "boost::math::find_scale<Dist, Policy>(%1%, %1%, %1%, Policy)"; 47 48 if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1)) 49 { 50 return policies::raise_domain_error<typename Dist::value_type>( 51 function, "Probability parameter was %1%, but must be >= 0 and <= 1!", p, pol); 52 } 53 if(!(boost::math::isfinite)(z)) 54 { 55 return policies::raise_domain_error<typename Dist::value_type>( 56 function, "find_scale z parameter was %1%, but must be finite!", z, pol); 57 } 58 if(!(boost::math::isfinite)(location)) 59 { 60 return policies::raise_domain_error<typename Dist::value_type>( 61 function, "find_scale location parameter was %1%, but must be finite!", location, pol); 62 } 63 64 //cout << "z " << z << ", p " << p << ", quantile(Dist(), p) " 65 //<< quantile(Dist(), p) << ", z - mean " << z - location 66 //<<", sd " << (z - location) / quantile(Dist(), p) << endl; 67 68 //quantile(N01, 0.001) -3.09023 69 //quantile(N01, 0.01) -2.32635 70 //quantile(N01, 0.05) -1.64485 71 //quantile(N01, 0.333333) -0.430728 72 //quantile(N01, 0.5) 0 73 //quantile(N01, 0.666667) 0.430728 74 //quantile(N01, 0.9) 1.28155 75 //quantile(N01, 0.95) 1.64485 76 //quantile(N01, 0.99) 2.32635 77 //quantile(N01, 0.999) 3.09023 78 79 typename Dist::value_type result = 80 (z - location) // difference between desired x and current location. 81 / quantile(Dist(), p); // standard distribution. 82 83 if (result <= 0) 84 { // If policy isn't to throw, return the scale <= 0. 85 policies::raise_evaluation_error<typename Dist::value_type>(function, 86 "Computed scale (%1%) is <= 0!" " Was the complement intended?", 87 result, Policy()); 88 } 89 return result; 90 } // template <class Dist, class Policy> find_scale 91 92 template <class Dist> 93 inline // with default policy. find_scale(typename Dist::value_type z,typename Dist::value_type p,typename Dist::value_type location)94 typename Dist::value_type find_scale( // For example, normal mean. 95 typename Dist::value_type z, // location of random variable z to give probability, P(X > z) == p. 96 // For example, a nominal minimum acceptable z, so that p * 100 % are > z 97 typename Dist::value_type p, // probability value desired at x, say 0.95 for 95% > z. 98 typename Dist::value_type location) // location parameter, for example, mean. 99 { // Forward to find_scale using the default policy. 100 return (find_scale<Dist>(z, p, location, policies::policy<>())); 101 } // find_scale 102 103 template <class Dist, class Real1, class Real2, class Real3, class Policy> find_scale(complemented4_type<Real1,Real2,Real3,Policy> const & c)104 inline typename Dist::value_type find_scale( 105 complemented4_type<Real1, Real2, Real3, Policy> const& c) 106 { 107 //cout << "cparam1 q " << c.param1 // q 108 // << ", c.dist z " << c.dist // z 109 // << ", c.param2 l " << c.param2 // l 110 // << ", quantile (Dist(), c.param1 = q) " 111 // << quantile(Dist(), c.param1) //q 112 // << endl; 113 114 #if !defined(BOOST_NO_SFINAE) && !BOOST_WORKAROUND(__SUNPRO_CC, BOOST_TESTED_AT(0x590)) 115 BOOST_STATIC_ASSERT(::boost::math::tools::is_distribution<Dist>::value); 116 BOOST_STATIC_ASSERT(::boost::math::tools::is_scaled_distribution<Dist>::value); 117 #endif 118 static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))"; 119 120 // Checks on arguments, as not complemented version, 121 // Explicit policy. 122 typename Dist::value_type q = c.param1; 123 if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1)) 124 { 125 return policies::raise_domain_error<typename Dist::value_type>( 126 function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, c.param3); 127 } 128 typename Dist::value_type z = c.dist; 129 if(!(boost::math::isfinite)(z)) 130 { 131 return policies::raise_domain_error<typename Dist::value_type>( 132 function, "find_scale z parameter was %1%, but must be finite!", z, c.param3); 133 } 134 typename Dist::value_type location = c.param2; 135 if(!(boost::math::isfinite)(location)) 136 { 137 return policies::raise_domain_error<typename Dist::value_type>( 138 function, "find_scale location parameter was %1%, but must be finite!", location, c.param3); 139 } 140 141 typename Dist::value_type result = 142 (c.dist - c.param2) // difference between desired x and current location. 143 / quantile(complement(Dist(), c.param1)); 144 // ( z - location) / (quantile(complement(Dist(), q)) 145 if (result <= 0) 146 { // If policy isn't to throw, return the scale <= 0. 147 policies::raise_evaluation_error<typename Dist::value_type>(function, 148 "Computed scale (%1%) is <= 0!" " Was the complement intended?", 149 result, Policy()); 150 } 151 return result; 152 } // template <class Dist, class Policy, class Real1, class Real2, class Real3> typename Dist::value_type find_scale 153 154 // So the user can start from the complement q = (1 - p) of the probability p, 155 // for example, s = find_scale<normal>(complement(z, q, l)); 156 157 template <class Dist, class Real1, class Real2, class Real3> find_scale(complemented3_type<Real1,Real2,Real3> const & c)158 inline typename Dist::value_type find_scale( 159 complemented3_type<Real1, Real2, Real3> const& c) 160 { 161 //cout << "cparam1 q " << c.param1 // q 162 // << ", c.dist z " << c.dist // z 163 // << ", c.param2 l " << c.param2 // l 164 // << ", quantile (Dist(), c.param1 = q) " 165 // << quantile(Dist(), c.param1) //q 166 // << endl; 167 168 #if !defined(BOOST_NO_SFINAE) && !BOOST_WORKAROUND(__SUNPRO_CC, BOOST_TESTED_AT(0x590)) 169 BOOST_STATIC_ASSERT(::boost::math::tools::is_distribution<Dist>::value); 170 BOOST_STATIC_ASSERT(::boost::math::tools::is_scaled_distribution<Dist>::value); 171 #endif 172 static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))"; 173 174 // Checks on arguments, as not complemented version, 175 // default policy policies::policy<>(). 176 typename Dist::value_type q = c.param1; 177 if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1)) 178 { 179 return policies::raise_domain_error<typename Dist::value_type>( 180 function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, policies::policy<>()); 181 } 182 typename Dist::value_type z = c.dist; 183 if(!(boost::math::isfinite)(z)) 184 { 185 return policies::raise_domain_error<typename Dist::value_type>( 186 function, "find_scale z parameter was %1%, but must be finite!", z, policies::policy<>()); 187 } 188 typename Dist::value_type location = c.param2; 189 if(!(boost::math::isfinite)(location)) 190 { 191 return policies::raise_domain_error<typename Dist::value_type>( 192 function, "find_scale location parameter was %1%, but must be finite!", location, policies::policy<>()); 193 } 194 195 typename Dist::value_type result = 196 (z - location) // difference between desired x and current location. 197 / quantile(complement(Dist(), q)); 198 // ( z - location) / (quantile(complement(Dist(), q)) 199 if (result <= 0) 200 { // If policy isn't to throw, return the scale <= 0. 201 policies::raise_evaluation_error<typename Dist::value_type>(function, 202 "Computed scale (%1%) is <= 0!" " Was the complement intended?", 203 result, policies::policy<>()); // This is only the default policy - also Want a version with Policy here. 204 } 205 return result; 206 } // template <class Dist, class Real1, class Real2, class Real3> typename Dist::value_type find_scale 207 208 } // namespace boost 209 } // namespace math 210 211 #endif // BOOST_STATS_FIND_SCALE_HPP 212