]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // find_location.cpp |
2 | ||
3 | // Copyright Paul A. Bristow 2008, 2010. | |
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 | // Example of finding location (mean) | |
11 | // for normal (Gaussian) & Cauchy distribution. | |
12 | ||
13 | // Note that this file contains Quickbook mark-up as well as code | |
14 | // and comments, don't change any of the special comment mark-ups! | |
15 | ||
16 | //#ifdef _MSC_VER | |
17 | //# pragma warning(disable: 4180) // qualifier has no effect (in Fusion). | |
18 | //#endif | |
19 | ||
20 | //[find_location1 | |
21 | /*` | |
22 | First we need some includes to access the normal distribution, | |
23 | the algorithms to find location (and some std output of course). | |
24 | */ | |
25 | ||
26 | #include <boost/math/distributions/normal.hpp> // for normal_distribution | |
27 | using boost::math::normal; // typedef provides default type is double. | |
28 | #include <boost/math/distributions/cauchy.hpp> // for cauchy_distribution | |
29 | using boost::math::cauchy; // typedef provides default type is double. | |
30 | #include <boost/math/distributions/find_location.hpp> | |
31 | using boost::math::find_location; // for mean | |
32 | #include <boost/math/distributions/find_scale.hpp> | |
33 | using boost::math::find_scale; // for standard devation | |
34 | using boost::math::complement; // Needed if you want to use the complement version. | |
35 | using boost::math::policies::policy; | |
36 | ||
37 | #include <iostream> | |
38 | using std::cout; using std::endl; | |
39 | #include <iomanip> | |
40 | using std::setw; using std::setprecision; | |
41 | #include <limits> | |
42 | using std::numeric_limits; | |
43 | ||
44 | //] [/find_location1] | |
45 | ||
46 | int main() | |
47 | { | |
48 | cout << "Example: Find location (or mean)." << endl; | |
49 | try | |
50 | { | |
51 | //[find_location2 | |
52 | /*` | |
53 | For this example, we will use the standard normal distribution, | |
54 | with mean (location) zero and standard deviation (scale) unity. | |
55 | This is also the default for this implementation. | |
56 | */ | |
57 | normal N01; // Default 'standard' normal distribution with zero mean and | |
58 | double sd = 1.; // normal default standard deviation is 1. | |
59 | /*`Suppose we want to find a different normal distribution whose mean is shifted | |
60 | so that only fraction p (here 0.001 or 0.1%) are below a certain chosen limit | |
61 | (here -2, two standard deviations). | |
62 | */ | |
63 | double z = -2.; // z to give prob p | |
64 | double p = 0.001; // only 0.1% below z | |
65 | ||
66 | cout << "Normal distribution with mean = " << N01.location() | |
67 | << ", standard deviation " << N01.scale() | |
68 | << ", has " << "fraction <= " << z | |
69 | << ", p = " << cdf(N01, z) << endl; | |
70 | cout << "Normal distribution with mean = " << N01.location() | |
71 | << ", standard deviation " << N01.scale() | |
72 | << ", has " << "fraction > " << z | |
73 | << ", p = " << cdf(complement(N01, z)) << endl; // Note: uses complement. | |
74 | /*` | |
75 | [pre | |
76 | Normal distribution with mean = 0, standard deviation 1, has fraction <= -2, p = 0.0227501 | |
77 | Normal distribution with mean = 0, standard deviation 1, has fraction > -2, p = 0.97725 | |
78 | ] | |
79 | We can now use ''find_location'' to give a new offset mean. | |
80 | */ | |
81 | double l = find_location<normal>(z, p, sd); | |
82 | cout << "offset location (mean) = " << l << endl; | |
83 | /*` | |
84 | that outputs: | |
85 | [pre | |
86 | offset location (mean) = 1.09023 | |
87 | ] | |
88 | showing that we need to shift the mean just over one standard deviation from its previous value of zero. | |
89 | ||
90 | Then we can check that we have achieved our objective | |
91 | by constructing a new distribution | |
92 | with the offset mean (but same standard deviation): | |
93 | */ | |
94 | normal np001pc(l, sd); // Same standard_deviation (scale) but with mean (location) shifted. | |
95 | /*` | |
96 | And re-calculating the fraction below our chosen limit. | |
97 | */ | |
98 | cout << "Normal distribution with mean = " << l | |
99 | << " has " << "fraction <= " << z | |
100 | << ", p = " << cdf(np001pc, z) << endl; | |
101 | cout << "Normal distribution with mean = " << l | |
102 | << " has " << "fraction > " << z | |
103 | << ", p = " << cdf(complement(np001pc, z)) << endl; | |
104 | /*` | |
105 | [pre | |
106 | Normal distribution with mean = 1.09023 has fraction <= -2, p = 0.001 | |
107 | Normal distribution with mean = 1.09023 has fraction > -2, p = 0.999 | |
108 | ] | |
109 | ||
110 | [h4 Controlling Error Handling from find_location] | |
111 | We can also control the policy for handling various errors. | |
112 | For example, we can define a new (possibly unwise) | |
113 | policy to ignore domain errors ('bad' arguments). | |
114 | ||
115 | Unless we are using the boost::math namespace, we will need: | |
116 | */ | |
117 | using boost::math::policies::policy; | |
118 | using boost::math::policies::domain_error; | |
119 | using boost::math::policies::ignore_error; | |
120 | ||
121 | /*` | |
122 | Using a typedef is often convenient, especially if it is re-used, | |
123 | although it is not required, as the various examples below show. | |
124 | */ | |
125 | typedef policy<domain_error<ignore_error> > ignore_domain_policy; | |
126 | // find_location with new policy, using typedef. | |
127 | l = find_location<normal>(z, p, sd, ignore_domain_policy()); | |
128 | // Default policy policy<>, needs "using boost::math::policies::policy;" | |
129 | l = find_location<normal>(z, p, sd, policy<>()); | |
130 | // Default policy, fully specified. | |
131 | l = find_location<normal>(z, p, sd, boost::math::policies::policy<>()); | |
132 | // A new policy, ignoring domain errors, without using a typedef. | |
133 | l = find_location<normal>(z, p, sd, policy<domain_error<ignore_error> >()); | |
134 | /*` | |
135 | If we want to use a probability that is the __complements of our probability, | |
136 | we should not even think of writing `find_location<normal>(z, 1 - p, sd)`, | |
137 | but use the complement version, see __why_complements. | |
138 | */ | |
139 | z = 2.; | |
140 | double q = 0.95; // = 1 - p; // complement. | |
141 | l = find_location<normal>(complement(z, q, sd)); | |
142 | ||
143 | normal np95pc(l, sd); // Same standard_deviation (scale) but with mean(location) shifted | |
144 | cout << "Normal distribution with mean = " << l << " has " | |
145 | << "fraction <= " << z << " = " << cdf(np95pc, z) << endl; | |
146 | cout << "Normal distribution with mean = " << l << " has " | |
147 | << "fraction > " << z << " = " << cdf(complement(np95pc, z)) << endl; | |
148 | //] [/find_location2] | |
149 | } | |
150 | catch(const std::exception& e) | |
151 | { // Always useful to include try & catch blocks because default policies | |
152 | // are to throw exceptions on arguments that cause errors like underflow, overflow. | |
153 | // Lacking try & catch blocks, the program will abort without a message below, | |
154 | // which may give some helpful clues as to the cause of the exception. | |
155 | std::cout << | |
156 | "\n""Message from thrown exception was:\n " << e.what() << std::endl; | |
157 | } | |
158 | return 0; | |
159 | } // int main() | |
160 | ||
161 | //[find_location_example_output | |
162 | /*` | |
163 | [pre | |
164 | Example: Find location (mean). | |
165 | Normal distribution with mean = 0, standard deviation 1, has fraction <= -2, p = 0.0227501 | |
166 | Normal distribution with mean = 0, standard deviation 1, has fraction > -2, p = 0.97725 | |
167 | offset location (mean) = 1.09023 | |
168 | Normal distribution with mean = 1.09023 has fraction <= -2, p = 0.001 | |
169 | Normal distribution with mean = 1.09023 has fraction > -2, p = 0.999 | |
170 | Normal distribution with mean = 0.355146 has fraction <= 2 = 0.95 | |
171 | Normal distribution with mean = 0.355146 has fraction > 2 = 0.05 | |
172 | ] | |
173 | */ | |
174 | //] [/find_location_example_output] |