]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | |
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. | |
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> | |
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> | |
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 |