find_scale.hpp 9.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211
  1. // Copyright John Maddock 2007.
  2. // Copyright Paul A. Bristow 2007.
  3. // Use, modification and distribution are subject to the
  4. // Boost Software License, Version 1.0. (See accompanying file
  5. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  6. #ifndef BOOST_STATS_FIND_SCALE_HPP
  7. #define BOOST_STATS_FIND_SCALE_HPP
  8. #include <boost/math/distributions/fwd.hpp> // for all distribution signatures.
  9. #include <boost/math/distributions/complement.hpp>
  10. #include <boost/math/policies/policy.hpp>
  11. // using boost::math::policies::policy;
  12. #include <boost/math/tools/traits.hpp>
  13. #include <boost/static_assert.hpp>
  14. #include <boost/math/special_functions/fpclassify.hpp>
  15. #include <boost/math/policies/error_handling.hpp>
  16. // using boost::math::complement; // will be needed by users who want complement,
  17. // but NOT placed here to avoid putting it in global scope.
  18. namespace boost
  19. {
  20. namespace math
  21. {
  22. // Function to find location of random variable z
  23. // to give probability p (given scale)
  24. // Applies to normal, lognormal, extreme value, Cauchy, (and symmetrical triangular),
  25. // distributions that have scale.
  26. // BOOST_STATIC_ASSERTs, see below, are used to enforce this.
  27. template <class Dist, class Policy>
  28. inline
  29. typename Dist::value_type find_scale( // For example, normal mean.
  30. typename Dist::value_type z, // location of random variable z to give probability, P(X > z) == p.
  31. // For example, a nominal minimum acceptable weight z, so that p * 100 % are > z
  32. typename Dist::value_type p, // probability value desired at x, say 0.95 for 95% > z.
  33. typename Dist::value_type location, // location parameter, for example, normal distribution mean.
  34. const Policy& pol
  35. )
  36. {
  37. #if !defined(BOOST_NO_SFINAE) && !BOOST_WORKAROUND(__SUNPRO_CC, BOOST_TESTED_AT(0x590))
  38. BOOST_STATIC_ASSERT(::boost::math::tools::is_distribution<Dist>::value);
  39. BOOST_STATIC_ASSERT(::boost::math::tools::is_scaled_distribution<Dist>::value);
  40. #endif
  41. static const char* function = "boost::math::find_scale<Dist, Policy>(%1%, %1%, %1%, Policy)";
  42. if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1))
  43. {
  44. return policies::raise_domain_error<typename Dist::value_type>(
  45. function, "Probability parameter was %1%, but must be >= 0 and <= 1!", p, pol);
  46. }
  47. if(!(boost::math::isfinite)(z))
  48. {
  49. return policies::raise_domain_error<typename Dist::value_type>(
  50. function, "find_scale z parameter was %1%, but must be finite!", z, pol);
  51. }
  52. if(!(boost::math::isfinite)(location))
  53. {
  54. return policies::raise_domain_error<typename Dist::value_type>(
  55. function, "find_scale location parameter was %1%, but must be finite!", location, pol);
  56. }
  57. //cout << "z " << z << ", p " << p << ", quantile(Dist(), p) "
  58. //<< quantile(Dist(), p) << ", z - mean " << z - location
  59. //<<", sd " << (z - location) / quantile(Dist(), p) << endl;
  60. //quantile(N01, 0.001) -3.09023
  61. //quantile(N01, 0.01) -2.32635
  62. //quantile(N01, 0.05) -1.64485
  63. //quantile(N01, 0.333333) -0.430728
  64. //quantile(N01, 0.5) 0
  65. //quantile(N01, 0.666667) 0.430728
  66. //quantile(N01, 0.9) 1.28155
  67. //quantile(N01, 0.95) 1.64485
  68. //quantile(N01, 0.99) 2.32635
  69. //quantile(N01, 0.999) 3.09023
  70. typename Dist::value_type result =
  71. (z - location) // difference between desired x and current location.
  72. / quantile(Dist(), p); // standard distribution.
  73. if (result <= 0)
  74. { // If policy isn't to throw, return the scale <= 0.
  75. policies::raise_evaluation_error<typename Dist::value_type>(function,
  76. "Computed scale (%1%) is <= 0!" " Was the complement intended?",
  77. result, Policy());
  78. }
  79. return result;
  80. } // template <class Dist, class Policy> find_scale
  81. template <class Dist>
  82. inline // with default policy.
  83. typename Dist::value_type find_scale( // For example, normal mean.
  84. typename Dist::value_type z, // location of random variable z to give probability, P(X > z) == p.
  85. // For example, a nominal minimum acceptable z, so that p * 100 % are > z
  86. typename Dist::value_type p, // probability value desired at x, say 0.95 for 95% > z.
  87. typename Dist::value_type location) // location parameter, for example, mean.
  88. { // Forward to find_scale using the default policy.
  89. return (find_scale<Dist>(z, p, location, policies::policy<>()));
  90. } // find_scale
  91. template <class Dist, class Real1, class Real2, class Real3, class Policy>
  92. inline typename Dist::value_type find_scale(
  93. complemented4_type<Real1, Real2, Real3, Policy> const& c)
  94. {
  95. //cout << "cparam1 q " << c.param1 // q
  96. // << ", c.dist z " << c.dist // z
  97. // << ", c.param2 l " << c.param2 // l
  98. // << ", quantile (Dist(), c.param1 = q) "
  99. // << quantile(Dist(), c.param1) //q
  100. // << endl;
  101. #if !defined(BOOST_NO_SFINAE) && !BOOST_WORKAROUND(__SUNPRO_CC, BOOST_TESTED_AT(0x590))
  102. BOOST_STATIC_ASSERT(::boost::math::tools::is_distribution<Dist>::value);
  103. BOOST_STATIC_ASSERT(::boost::math::tools::is_scaled_distribution<Dist>::value);
  104. #endif
  105. static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))";
  106. // Checks on arguments, as not complemented version,
  107. // Explicit policy.
  108. typename Dist::value_type q = c.param1;
  109. if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1))
  110. {
  111. return policies::raise_domain_error<typename Dist::value_type>(
  112. function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, c.param3);
  113. }
  114. typename Dist::value_type z = c.dist;
  115. if(!(boost::math::isfinite)(z))
  116. {
  117. return policies::raise_domain_error<typename Dist::value_type>(
  118. function, "find_scale z parameter was %1%, but must be finite!", z, c.param3);
  119. }
  120. typename Dist::value_type location = c.param2;
  121. if(!(boost::math::isfinite)(location))
  122. {
  123. return policies::raise_domain_error<typename Dist::value_type>(
  124. function, "find_scale location parameter was %1%, but must be finite!", location, c.param3);
  125. }
  126. typename Dist::value_type result =
  127. (c.dist - c.param2) // difference between desired x and current location.
  128. / quantile(complement(Dist(), c.param1));
  129. // ( z - location) / (quantile(complement(Dist(), q))
  130. if (result <= 0)
  131. { // If policy isn't to throw, return the scale <= 0.
  132. policies::raise_evaluation_error<typename Dist::value_type>(function,
  133. "Computed scale (%1%) is <= 0!" " Was the complement intended?",
  134. result, Policy());
  135. }
  136. return result;
  137. } // template <class Dist, class Policy, class Real1, class Real2, class Real3> typename Dist::value_type find_scale
  138. // So the user can start from the complement q = (1 - p) of the probability p,
  139. // for example, s = find_scale<normal>(complement(z, q, l));
  140. template <class Dist, class Real1, class Real2, class Real3>
  141. inline typename Dist::value_type find_scale(
  142. complemented3_type<Real1, Real2, Real3> const& c)
  143. {
  144. //cout << "cparam1 q " << c.param1 // q
  145. // << ", c.dist z " << c.dist // z
  146. // << ", c.param2 l " << c.param2 // l
  147. // << ", quantile (Dist(), c.param1 = q) "
  148. // << quantile(Dist(), c.param1) //q
  149. // << endl;
  150. #if !defined(BOOST_NO_SFINAE) && !BOOST_WORKAROUND(__SUNPRO_CC, BOOST_TESTED_AT(0x590))
  151. BOOST_STATIC_ASSERT(::boost::math::tools::is_distribution<Dist>::value);
  152. BOOST_STATIC_ASSERT(::boost::math::tools::is_scaled_distribution<Dist>::value);
  153. #endif
  154. static const char* function = "boost::math::find_scale<Dist, Policy>(complement(%1%, %1%, %1%, Policy))";
  155. // Checks on arguments, as not complemented version,
  156. // default policy policies::policy<>().
  157. typename Dist::value_type q = c.param1;
  158. if(!(boost::math::isfinite)(q) || (q < 0) || (q > 1))
  159. {
  160. return policies::raise_domain_error<typename Dist::value_type>(
  161. function, "Probability parameter was %1%, but must be >= 0 and <= 1!", q, policies::policy<>());
  162. }
  163. typename Dist::value_type z = c.dist;
  164. if(!(boost::math::isfinite)(z))
  165. {
  166. return policies::raise_domain_error<typename Dist::value_type>(
  167. function, "find_scale z parameter was %1%, but must be finite!", z, policies::policy<>());
  168. }
  169. typename Dist::value_type location = c.param2;
  170. if(!(boost::math::isfinite)(location))
  171. {
  172. return policies::raise_domain_error<typename Dist::value_type>(
  173. function, "find_scale location parameter was %1%, but must be finite!", location, policies::policy<>());
  174. }
  175. typename Dist::value_type result =
  176. (z - location) // difference between desired x and current location.
  177. / quantile(complement(Dist(), q));
  178. // ( z - location) / (quantile(complement(Dist(), q))
  179. if (result <= 0)
  180. { // If policy isn't to throw, return the scale <= 0.
  181. policies::raise_evaluation_error<typename Dist::value_type>(function,
  182. "Computed scale (%1%) is <= 0!" " Was the complement intended?",
  183. result, policies::policy<>()); // This is only the default policy - also Want a version with Policy here.
  184. }
  185. return result;
  186. } // template <class Dist, class Real1, class Real2, class Real3> typename Dist::value_type find_scale
  187. } // namespace boost
  188. } // namespace math
  189. #endif // BOOST_STATS_FIND_SCALE_HPP