]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // distribution_construction.cpp |
2 | ||
3 | // Copyright Paul A. Bristow 2007, 2010, 2012. | |
4 | ||
5 | // Use, modification and distribution are subject to the | |
6 | // Boost Software License, Version 1.0. | |
7 | // (See accompanying file LICENSE_1_0.txt | |
8 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
9 | ||
10 | // Caution: this file contains Quickbook markup as well as code | |
11 | // and comments, don't change any of the special comment markups! | |
12 | ||
13 | #ifdef _MSC_VER | |
14 | # pragma warning (disable : 4996) // disable -D_SCL_SECURE_NO_WARNINGS C++ 'Checked Iterators' | |
15 | #endif | |
16 | ||
17 | #include <iostream> | |
18 | #include <exception> | |
19 | ||
20 | //[distribution_construction_1 | |
21 | ||
22 | /*` | |
23 | The structure of distributions is rather different from some other statistical libraries, | |
92f5a8d4 TL |
24 | for example, those written in less object-oriented language like FORTRAN and C that |
25 | provide a few arguments to each free function. | |
7c673cae | 26 | |
92f5a8d4 | 27 | Boost.Math library instead provides each distribution as a template C++ class. |
7c673cae FG |
28 | A distribution is constructed with a few arguments, and then |
29 | member and non-member functions are used to find values of the | |
30 | distribution, often a function of a random variate. | |
31 | ||
32 | For this demonstration, first we need some includes to access the | |
33 | negative binomial distribution (and the binomial, beta and gamma distributions too). | |
34 | ||
35 | To demonstrate the use with a high precision User-defined floating-point type | |
92f5a8d4 TL |
36 | `cpp_bin_float`, we also need an include from Boost.Multiprecision. |
37 | (We could equally well have used a `cpp_dec_float` multiprecision type). | |
38 | ||
39 | We choose a typedef `cpp_bin_float_50` to provide a 50 decimal digit type, | |
40 | but we could equally have chosen at 128-bit type `cpp_bin_float_quad`, | |
41 | or on some platforms `__float128`, providing about 35 decimal digits. | |
7c673cae FG |
42 | */ |
43 | ||
44 | #include <boost/math/distributions/negative_binomial.hpp> // for negative_binomial_distribution | |
45 | using boost::math::negative_binomial_distribution; // default type is double. | |
46 | using boost::math::negative_binomial; // typedef provides default type is double. | |
47 | #include <boost/math/distributions/binomial.hpp> // for binomial_distribution. | |
48 | #include <boost/math/distributions/beta.hpp> // for beta_distribution. | |
49 | #include <boost/math/distributions/gamma.hpp> // for gamma_distribution. | |
50 | #include <boost/math/distributions/normal.hpp> // for normal_distribution. | |
51 | ||
92f5a8d4 | 52 | #include <boost/multiprecision/cpp_bin_float.hpp> // for cpp_bin_float_50 |
7c673cae FG |
53 | /*` |
54 | Several examples of constructing distributions follow: | |
55 | */ | |
56 | //] [/distribution_construction_1 end of Quickbook in C++ markup] | |
57 | ||
58 | int main() | |
59 | { | |
60 | try | |
61 | { | |
62 | //[distribution_construction_2 | |
63 | /*` | |
64 | First, a negative binomial distribution with 8 successes | |
65 | and a success fraction 0.25, 25% or 1 in 4, is constructed like this: | |
66 | */ | |
67 | boost::math::negative_binomial_distribution<double> mydist0(8., 0.25); | |
68 | /*` | |
69 | But this is inconveniently long, so we might be tempted to write | |
70 | */ | |
71 | using namespace boost::math; | |
72 | /*` | |
73 | but this might risk ambiguity with names in `std random` so | |
74 | [*much] better is explicit `using boost::math::` statements, for example: | |
75 | */ | |
76 | using boost::math::negative_binomial_distribution; | |
77 | /*` | |
78 | and we can still reduce typing. | |
79 | ||
80 | Since the vast majority of applications use will be using `double` precision, | |
81 | the template argument to the distribution (`RealType`) defaults | |
82 | to type `double`, so we can also write: | |
83 | */ | |
84 | ||
85 | negative_binomial_distribution<> mydist9(8., 0.25); // Uses default `RealType = double`. | |
86 | ||
87 | /*` | |
88 | But the name `negative_binomial_distribution` is still inconveniently long, | |
89 | so, for most distributions, a convenience `typedef` is provided, for example: | |
90 | ||
91 | typedef negative_binomial_distribution<double> negative_binomial; // Reserved name of type double. | |
92 | ||
93 | [caution | |
94 | This convenience typedef is [*not provided] if a clash would occur | |
92f5a8d4 | 95 | with the name of a function; currently only `beta` and `gamma` |
7c673cae FG |
96 | fall into this category. |
97 | ] | |
98 | ||
99 | So, after a using statement, | |
100 | */ | |
101 | ||
102 | using boost::math::negative_binomial; | |
103 | ||
104 | /*` | |
105 | we have a convenient typedef to `negative_binomial_distribution<double>`: | |
106 | */ | |
107 | negative_binomial mydist(8., 0.25); | |
108 | ||
109 | /*` | |
110 | Some more examples using the convenience typedef: | |
111 | */ | |
112 | negative_binomial mydist10(5., 0.4); // Both arguments double. | |
113 | /*` | |
92f5a8d4 | 114 | And automatic conversion of arguments takes place, so you can use integers and floats: |
7c673cae | 115 | */ |
92f5a8d4 | 116 | negative_binomial mydist11(5, 0.4); // Using provided typedef of type double, and int and double arguments. |
7c673cae FG |
117 | /*` |
118 | This is probably the most common usage. | |
92f5a8d4 | 119 | Other combination are possible too: |
7c673cae FG |
120 | */ |
121 | negative_binomial mydist12(5., 0.4F); // Double and float arguments. | |
122 | negative_binomial mydist13(5, 1); // Both arguments integer. | |
123 | ||
124 | /*` | |
125 | Similarly for most other distributions like the binomial. | |
126 | */ | |
127 | binomial mybinomial(1, 0.5); // is more concise than | |
128 | binomial_distribution<> mybinomd1(1, 0.5); | |
129 | ||
130 | /*` | |
131 | For cases when the typdef distribution name would clash with a math special function | |
132 | (currently only beta and gamma) | |
133 | the typedef is deliberately not provided, and the longer version of the name | |
92f5a8d4 | 134 | must be used, so for example, do not use: |
7c673cae FG |
135 | |
136 | using boost::math::beta; | |
137 | beta mybetad0(1, 0.5); // Error beta is a math FUNCTION! | |
138 | ||
139 | Which produces the error messages: | |
140 | ||
141 | [pre | |
142 | error C2146: syntax error : missing ';' before identifier 'mybetad0' | |
143 | warning C4551: function call missing argument list | |
144 | error C3861: 'mybetad0': identifier not found | |
145 | ] | |
146 | ||
147 | Instead you should use: | |
148 | */ | |
149 | using boost::math::beta_distribution; | |
150 | beta_distribution<> mybetad1(1, 0.5); | |
151 | /*` | |
152 | or for the gamma distribution: | |
153 | */ | |
154 | gamma_distribution<> mygammad1(1, 0.5); | |
155 | ||
156 | /*` | |
157 | We can, of course, still provide the type explicitly thus: | |
158 | */ | |
159 | ||
160 | // Explicit double precision: both arguments are double: | |
161 | negative_binomial_distribution<double> mydist1(8., 0.25); | |
162 | ||
163 | // Explicit float precision, double arguments are truncated to float: | |
164 | negative_binomial_distribution<float> mydist2(8., 0.25); | |
165 | ||
166 | // Explicit float precision, integer & double arguments converted to float: | |
167 | negative_binomial_distribution<float> mydist3(8, 0.25); | |
168 | ||
169 | // Explicit float precision, float arguments, so no conversion: | |
170 | negative_binomial_distribution<float> mydist4(8.F, 0.25F); | |
171 | ||
172 | // Explicit float precision, integer arguments promoted to float. | |
173 | negative_binomial_distribution<float> mydist5(8, 1); | |
174 | ||
175 | // Explicit double precision: | |
92f5a8d4 | 176 | negative_binomial_distribution<double> mydist6(5., 0.4); |
7c673cae FG |
177 | |
178 | // Explicit long double precision: | |
179 | negative_binomial_distribution<long double> mydist7(8., 0.25); | |
92f5a8d4 | 180 | |
7c673cae | 181 | /*` |
92f5a8d4 TL |
182 | And you can use your own template RealType, |
183 | for example, `boost::math::cpp_bin_float_50` (an arbitrary 50 decimal digits precision type), | |
7c673cae FG |
184 | then we can write: |
185 | */ | |
186 | using namespace boost::multiprecision; | |
92f5a8d4 | 187 | negative_binomial_distribution<cpp_bin_float_50> mydist8(8, 0.25); |
7c673cae | 188 | |
7c673cae FG |
189 | // `integer` arguments are promoted to your RealType exactly, but |
190 | // `double` argument are converted to RealType, | |
92f5a8d4 TL |
191 | // most likely losing precision! |
192 | ||
193 | // So DON'T be tempted to write the 'obvious': | |
194 | negative_binomial_distribution<cpp_bin_float_50> mydist20(8, 0.23456789012345678901234567890); | |
195 | // to avoid truncation of second parameter to `0.2345678901234567` and loss of precision. | |
7c673cae | 196 | |
92f5a8d4 TL |
197 | // Instead pass a quoted decimal digit string: |
198 | negative_binomial_distribution<cpp_bin_float_50> mydist21(8, cpp_bin_float_50("0.23456789012345678901234567890") ); | |
7c673cae FG |
199 | |
200 | // Ensure that all potentially significant digits are shown. | |
92f5a8d4 TL |
201 | std::cout.precision(std::numeric_limits<cpp_bin_float_50>::digits10); |
202 | // | |
203 | cpp_bin_float_50 x("1.23456789012345678901234567890"); | |
7c673cae FG |
204 | std::cout << pdf(mydist8, x) << std::endl; |
205 | /*` showing 0.00012630010495970320103876754721976419438231705359935 | |
92f5a8d4 | 206 | 0.00012630010495970320103876754721976419438231528547467 |
7c673cae FG |
207 | |
208 | [warning When using multiprecision, it is all too easy to get accidental truncation!] | |
209 | ||
210 | For example, if you write | |
211 | */ | |
212 | std::cout << pdf(mydist8, 1.23456789012345678901234567890) << std::endl; | |
213 | /*` | |
214 | showing 0.00012630010495970318465064569310967179576805651692929, | |
215 | which is wrong at about the 17th decimal digit! | |
216 | ||
217 | This is because the value provided is truncated to a `double`, effectively | |
218 | `double x = 1.23456789012345678901234567890;` | |
219 | ||
220 | Then the now `double x` is passed to function `pdf`, | |
92f5a8d4 | 221 | and this truncated `double` value is finally promoted to `cpp_bin_float_50`. |
7c673cae FG |
222 | |
223 | Another way of quietly getting the wrong answer is to write: | |
224 | */ | |
92f5a8d4 | 225 | std::cout << pdf(mydist8, cpp_bin_float_50(1.23456789012345678901234567890)) << std::endl; |
7c673cae FG |
226 | /*` |
227 | A correct way from a multi-digit string value is | |
228 | */ | |
92f5a8d4 | 229 | std::cout << pdf(mydist8, cpp_bin_float_50("1.23456789012345678901234567890")) << std::endl; |
7c673cae FG |
230 | /*` |
231 | ||
232 | [tip Getting about 17 decimal digits followed by many zeros is often a sign of accidental truncation.] | |
233 | */ | |
234 | ||
235 | /*` | |
236 | [h4 Default arguments to distribution constructors.] | |
237 | ||
238 | Note that default constructor arguments are only provided for some distributions. | |
239 | So if you wrongly assume a default argument, you will get an error message, for example: | |
240 | ||
241 | negative_binomial_distribution<> mydist8; | |
242 | ||
243 | [pre error C2512 no appropriate default constructor available.] | |
244 | ||
245 | No default constructors are provided for the `negative binomial` distribution, | |
246 | because it is difficult to chose any sensible default values for this distribution. | |
247 | ||
248 | For other distributions, like the normal distribution, | |
249 | it is obviously very useful to provide 'standard' | |
250 | defaults for the mean (zero) and standard deviation (unity) thus: | |
251 | ||
252 | normal_distribution(RealType mean = 0, RealType sd = 1); | |
253 | ||
92f5a8d4 | 254 | So in this case we can more tersely write: |
7c673cae FG |
255 | */ |
256 | using boost::math::normal; | |
257 | ||
92f5a8d4 | 258 | normal norm1; // Standard normal distribution N[0,1]. |
7c673cae FG |
259 | normal norm2(2); // Mean = 2, std deviation = 1. |
260 | normal norm3(2, 3); // Mean = 2, std deviation = 3. | |
261 | ||
262 | } | |
263 | catch(std::exception &ex) | |
264 | { | |
265 | std::cout << ex.what() << std::endl; | |
266 | } | |
267 | ||
268 | return 0; | |
269 | } // int main() | |
270 | ||
271 | /*`There is no useful output from this demonstration program, of course. */ | |
272 | ||
273 | //] [/end of distribution_construction_2] | |
274 | ||
275 | /* | |
276 | //[distribution_construction_output | |
277 | ||
278 | 0.00012630010495970320103876754721976419438231705359935 | |
279 | 0.00012630010495970318465064569310967179576805651692929 | |
280 | 0.00012630010495970318465064569310967179576805651692929 | |
281 | 0.00012630010495970320103876754721976419438231705359935 | |
282 | ||
283 | //] [/distribution_construction_output] | |
284 | ||
92f5a8d4 TL |
285 | |
286 | 0.00012630010495970320103876754721976419438231528547467 | |
287 | 0.0001263001049597031846506456931096717957680547488046 | |
288 | 0.0001263001049597031846506456931096717957680547488046 | |
289 | 0.00012630010495970320103876754721976419438231528547467 | |
290 | ||
291 | ||
7c673cae FG |
292 | */ |
293 | ||
294 | ||
295 |