neg_binom_confidence_limits.cpp 6.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178
  1. // neg_binomial_confidence_limits.cpp
  2. // Copyright John Maddock 2006
  3. // Copyright Paul A. Bristow 2007, 2010
  4. // Use, modification and distribution are subject to the
  5. // Boost Software License, Version 1.0.
  6. // (See accompanying file LICENSE_1_0.txt
  7. // or copy at http://www.boost.org/LICENSE_1_0.txt)
  8. // Caution: this file contains quickbook markup as well as code
  9. // and comments, don't change any of the special comment markups!
  10. //[neg_binomial_confidence_limits
  11. /*`
  12. First we need some includes to access the negative binomial distribution
  13. (and some basic std output of course).
  14. */
  15. #include <boost/math/distributions/negative_binomial.hpp>
  16. using boost::math::negative_binomial;
  17. #include <iostream>
  18. using std::cout; using std::endl;
  19. #include <iomanip>
  20. using std::setprecision;
  21. using std::setw; using std::left; using std::fixed; using std::right;
  22. /*`
  23. First define a table of significance levels: these are the
  24. probabilities that the true occurrence frequency lies outside the calculated
  25. interval:
  26. */
  27. double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
  28. /*`
  29. Confidence value as % is (1 - alpha) * 100, so alpha 0.05 == 95% confidence
  30. that the true occurrence frequency lies *inside* the calculated interval.
  31. We need a function to calculate and print confidence limits
  32. for an observed frequency of occurrence
  33. that follows a negative binomial distribution.
  34. */
  35. void confidence_limits_on_frequency(unsigned trials, unsigned successes)
  36. {
  37. // trials = Total number of trials.
  38. // successes = Total number of observed successes.
  39. // failures = trials - successes.
  40. // success_fraction = successes /trials.
  41. // Print out general info:
  42. cout <<
  43. "______________________________________________\n"
  44. "2-Sided Confidence Limits For Success Fraction\n"
  45. "______________________________________________\n\n";
  46. cout << setprecision(7);
  47. cout << setw(40) << left << "Number of trials" << " = " << trials << "\n";
  48. cout << setw(40) << left << "Number of successes" << " = " << successes << "\n";
  49. cout << setw(40) << left << "Number of failures" << " = " << trials - successes << "\n";
  50. cout << setw(40) << left << "Observed frequency of occurrence" << " = " << double(successes) / trials << "\n";
  51. // Print table header:
  52. cout << "\n\n"
  53. "___________________________________________\n"
  54. "Confidence Lower Upper\n"
  55. " Value (%) Limit Limit\n"
  56. "___________________________________________\n";
  57. /*`
  58. And now for the important part - the bounds themselves.
  59. For each value of /alpha/, we call `find_lower_bound_on_p` and
  60. `find_upper_bound_on_p` to obtain lower and upper bounds respectively.
  61. Note that since we are calculating a two-sided interval,
  62. we must divide the value of alpha in two. Had we been calculating a
  63. single-sided interval, for example: ['"Calculate a lower bound so that we are P%
  64. sure that the true occurrence frequency is greater than some value"]
  65. then we would *not* have divided by two.
  66. */
  67. // Now print out the upper and lower limits for the alpha table values.
  68. for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
  69. {
  70. // Confidence value:
  71. cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]);
  72. // Calculate bounds:
  73. double lower = negative_binomial::find_lower_bound_on_p(trials, successes, alpha[i]/2);
  74. double upper = negative_binomial::find_upper_bound_on_p(trials, successes, alpha[i]/2);
  75. // Print limits:
  76. cout << fixed << setprecision(5) << setw(15) << right << lower;
  77. cout << fixed << setprecision(5) << setw(15) << right << upper << endl;
  78. }
  79. cout << endl;
  80. } // void confidence_limits_on_frequency(unsigned trials, unsigned successes)
  81. /*`
  82. And then call confidence_limits_on_frequency with increasing numbers of trials,
  83. but always the same success fraction 0.1, or 1 in 10.
  84. */
  85. int main()
  86. {
  87. confidence_limits_on_frequency(20, 2); // 20 trials, 2 successes, 2 in 20, = 1 in 10 = 0.1 success fraction.
  88. confidence_limits_on_frequency(200, 20); // More trials, but same 0.1 success fraction.
  89. confidence_limits_on_frequency(2000, 200); // Many more trials, but same 0.1 success fraction.
  90. return 0;
  91. } // int main()
  92. //] [/negative_binomial_confidence_limits_eg end of Quickbook in C++ markup]
  93. /*
  94. ______________________________________________
  95. 2-Sided Confidence Limits For Success Fraction
  96. ______________________________________________
  97. Number of trials = 20
  98. Number of successes = 2
  99. Number of failures = 18
  100. Observed frequency of occurrence = 0.1
  101. ___________________________________________
  102. Confidence Lower Upper
  103. Value (%) Limit Limit
  104. ___________________________________________
  105. 50.000 0.04812 0.13554
  106. 75.000 0.03078 0.17727
  107. 90.000 0.01807 0.22637
  108. 95.000 0.01235 0.26028
  109. 99.000 0.00530 0.33111
  110. 99.900 0.00164 0.41802
  111. 99.990 0.00051 0.49202
  112. 99.999 0.00016 0.55574
  113. ______________________________________________
  114. 2-Sided Confidence Limits For Success Fraction
  115. ______________________________________________
  116. Number of trials = 200
  117. Number of successes = 20
  118. Number of failures = 180
  119. Observed frequency of occurrence = 0.1000000
  120. ___________________________________________
  121. Confidence Lower Upper
  122. Value (%) Limit Limit
  123. ___________________________________________
  124. 50.000 0.08462 0.11350
  125. 75.000 0.07580 0.12469
  126. 90.000 0.06726 0.13695
  127. 95.000 0.06216 0.14508
  128. 99.000 0.05293 0.16170
  129. 99.900 0.04343 0.18212
  130. 99.990 0.03641 0.20017
  131. 99.999 0.03095 0.21664
  132. ______________________________________________
  133. 2-Sided Confidence Limits For Success Fraction
  134. ______________________________________________
  135. Number of trials = 2000
  136. Number of successes = 200
  137. Number of failures = 1800
  138. Observed frequency of occurrence = 0.1000000
  139. ___________________________________________
  140. Confidence Lower Upper
  141. Value (%) Limit Limit
  142. ___________________________________________
  143. 50.000 0.09536 0.10445
  144. 75.000 0.09228 0.10776
  145. 90.000 0.08916 0.11125
  146. 95.000 0.08720 0.11352
  147. 99.000 0.08344 0.11802
  148. 99.900 0.07921 0.12336
  149. 99.990 0.07577 0.12795
  150. 99.999 0.07282 0.13206
  151. */