find_location_example.cpp 6.6 KB

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