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 deviation
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
main()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]
175