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