find_scale.hpp 9.9 KB

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