[Boost-commit] svn:boost r66721 - trunk/libs/math/example

Subject: [Boost-commit] svn:boost r66721 - trunk/libs/math/example
From: pbristow_at_[hidden]
Date: 2010-11-24 11:34:36

Author: pbristow
Date: 2010-11-24 11:34:31 EST (Wed, 24 Nov 2010)
New Revision: 66721
URL: http://svn.boost.org/trac/boost/changeset/66721

Log:
Added new distributions.
Added:
   trunk/libs/math/example/geometric_examples.cpp   (contents, props changed)
   trunk/libs/math/example/inverse_chi_squared_example.cpp   (contents, props changed)
   trunk/libs/math/example/inverse_chi_squared_find_df_example.cpp   (contents, props changed)
   trunk/libs/math/example/inverse_gamma_distribution_example.cpp   (contents, props changed)
   trunk/libs/math/example/inverse_gamma_example.cpp   (contents, props changed)

Added: trunk/libs/math/example/geometric_examples.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/math/example/geometric_examples.cpp	2010-11-24 11:34:31 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,354 @@
+// geometric_examples.cpp
+
+// Copyright Paul A. Bristow 2010.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// This file is written to be included from a Quickbook .qbk document.
+// It can still be compiled by the C++ compiler, and run. 
+// Any output can also be added here as comment or included or pasted in elsewhere.
+// Caution: this file contains Quickbook markup as well as code
+// and comments: don't change any of the special comment markups!
+
+// Examples of using the geometric distribution.
+
+//[geometric_eg1_1
+/*`
+For this example, we will opt to #define two macros to control 
+the error and discrete handling policies.
+For this simple example, we want to avoid throwing
+an exception (the default policy) and just return infinity.
+We want to treat the distribution as if it was continuous,
+so we choose a discrete_quantile policy of real,
+rather than the default policy integer_round_outwards.
+*/
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
+/*`
+[caution It is vital to #include distributions etc *after* the above #defines]
+After that we need some includes to provide easy access to the negative binomial distribution,
+and we need some std library iostream, of course.
+*/
+#include <boost/math/distributions/geometric.hpp>
+  // for geometric_distribution
+  using ::boost::math::geometric_distribution; // 
+  using ::boost::math::geometric; // typedef provides default type is double.
+  using  ::boost::math::pdf; // Probability mass function.
+  using  ::boost::math::cdf; // Cumulative density function.
+  using  ::boost::math::quantile;
+
+#include <boost/math/distributions/negative_binomial.hpp>
+  // for negative_binomial_distribution
+  using boost::math::negative_binomial; // typedef provides default type is double.
+
+#include <boost/math/distributions/normal.hpp>
+  // for negative_binomial_distribution
+  using boost::math::normal; // typedef provides default type is double.
+
+#include <iostream>
+  using std::cout; using std::endl;
+  using std::noshowpoint; using std::fixed; using std::right; using std::left;
+#include <iomanip>
+  using std::setprecision; using std::setw; 
+
+#include <limits>
+  using std::numeric_limits;
+//] [geometric_eg1_1]
+
+int main()
+{
+  cout <<"Geometric distribution example" << endl;
+  cout << endl;
+  cout.precision(std::numeric_limits<double>::max_digits10);
+  // gives all potentially significant digits for the type, here double.
+  cout.precision(4); 
+  try
+  {
+//[geometric_eg1_2
+/*`
+It is always sensible to use try and catch blocks because defaults policies are to
+throw an exception if anything goes wrong.
+
+Simple try'n'catch blocks (see below) will ensure that you get a
+helpful error message instead of an abrupt (and silent) program abort.
+
+[h6 Throwing a dice]
+The Geometric distribution describes the probability (/p/) of a number of failures
+to get the first success in /k/ Bernoulli trials.
+(A [@http://en.wikipedia.org/wiki/Bernoulli_distribution Bernoulli trial]
+is one with only two possible outcomes, success of failure,
+and /p/ is the probability of success).
+
+Suppose an 'fair' 6-face dice is thrown repeatedly: 
+*/
+    double success_fraction = 1./6; // success_fraction (p) = 0.1666
+    // (so failure_fraction is 1 - success_fraction = 5./6 = 1- 0.1666 = 0.8333)
+
+/*`If the dice is thrown repeatedly until the *first* time a /three/ appears.
+The probablility distribution of the number of times it is thrown *not* getting a /three/
+ (/not-a-threes/ number of failures to get a /three/)
+is a geometric distribution with the success_fraction = 1/6 = 0.1666[recur].
+
+We therefore start by constructing a geometric distribution
+with the one parameter success_fraction, the probability of success.
+*/
+    geometric g6(success_fraction); // type double by default.
+/*`
+To confirm, we can echo the success_fraction parameter of the distribution.
+*/
+    cout << "success fraction of a six-sided dice is " << g6.success_fraction() << endl;
+/*`So the probability of getting a three at the first throw (zero failures) is
+*/
+    cout << pdf(g6, 0) << endl; // 0.1667
+    cout << cdf(g6, 0) << endl; // 0.1667
+/*`Note that the cdf and pdf are identical because the is only one throw.
+If we want the probability of getting the first /three/ on the 2nd throw:
+*/
+    cout << pdf(g6, 1) << endl; // 0.1389
+
+/*`If we want the probability of getting the first /three/ on the 1st or 2nd throw
+(allowing one failure):
+*/
+    cout << "pdf(g6, 0) + pdf(g6, 1) = " << pdf(g6, 0) + pdf(g6, 1) << endl;
+/*`Or more conveniently, and more generally,
+we can use the Cumulative Distribution Function CDF.*/
+
+    cout << "cdf(g6, 1) = " << cdf(g6, 1) << endl; // 0.3056
+
+/*`If we allow many more (12) throws, the probability of getting our /three/ gets very high:*/
+    cout << "cdf(g6, 12) = " << cdf(g6, 12) << endl; // 0.9065 or 90% probability.
+/*`If we want to be much more confident, say 99%, 
+we can estimate the number of throws to be this sure
+using the inverse or quantile.
+*/
+    cout << "quantile(g6, 0.99) = " << quantile(g6, 0.99) << endl; // 24.26
+/*`Note that the value returned is not an integer:
+if you want an integer result you should use either floor, round or ceil functions,
+or use the policies mechanism.
+See [link math_toolkit.policy.pol_tutorial.understand_dis_quant
+Understanding Quantiles of Discrete Distributions] 
+
+The geometric distribution is related to the negative binomial
+__spaces `negative_binomial_distribution(RealType r, RealType p);` with parameter /r/ = 1.
+So we could get the same result using the negative binomial,
+but using the geometric the results will be faster, and may be more accurate.
+*/
+    negative_binomial nb(1, success_fraction);
+    cout << pdf(nb, 1) << endl; // 0.1389
+    cout << cdf(nb, 1) << endl; // 0.3056
+/*`We could also the complement to express the required probability
+as 1 - 0.99 = 0.01 (and get the same result):
+*/
+    cout << "quantile(complement(g6, 1 - p))  " << quantile(complement(g6, 0.01)) << endl; // 24.26
+/*`
+Note too that Boost.Math geometric distribution is implemented as a continuous function.
+Unlike other implementations (for example R) it *uses* the number of failures as a *real* parameter,
+not as an integer. If you want this integer behaviour, you may need to enforce this by
+rounding the parameter you pass, probably rounding down, to the nearest integer.
+For example, R returns the success fraction probability for all values of failures
+from 0 to 0.999999 thus:
+[pre
+__spaces R> formatC(pgeom(0.0001,0.5, FALSE), digits=17) "               0.5"
+] [/pre]
+So in Boost.Math the equivalent is
+*/
+    geometric g05(0.5);  // Probability of success = 0.5 or 50%
+    cout.precision(numeric_limits<double>::max_digits10); // Display all potentially significant digits.
+    cout << cdf(g05, 0.0001) << endl; // returns 0.5000346561579232, not exact 0.5.
+/*`To get the R discrete behaviour, you simply need to round with,
+for example, the `floor` function.
+*/
+    cout << cdf(g05, floor(0.0001)) << endl; // returns exactly 0.5
+/*`
+[pre
+`> formatC(pgeom(0.9999999,0.5, FALSE), digits=17) [1] "              0.25"`
+`> formatC(pgeom(1.999999,0.5, FALSE), digits=17)[1] "              0.25" k = 1`
+`> formatC(pgeom(1.9999999,0.5, FALSE), digits=17)[1] "0.12500000000000003" k = 2`
+] [/pre]
+shows that R makes an arbitrary round-up decision at about 1e7 from the next integer above.
+This may be convenient in practice, and could be replicated in C++ if desired.
+
+[h6 Surveying customers to find one with a faulty product]
+A company knows from warranty claims that 2% of their products will be faulty,
+so the 'success_fraction' of finding a fault is 0.02.
+It wants to interview a purchaser of faulty products to assess their 'user experience'.
+
+To estimate how many customers they will probably need to contact 
+in order to find one who has suffered from the fault,
+we first construct a geometric distribution with probability 0.02,
+and then chose a confidence, say 80%, 95%, or 99% to finding a customer with a fault.
+Finally, we probably want to round up the result to the integer above using the `ceil` function.
+(We could also use a policy, but that is hardly worthwhile for this simple application.)
+
+(This also assumes that each customer only buys one product:
+if customers bought more than one item,
+the probability of finding a customer with a fault obviously improves.)
+*/
+    cout.precision(5);
+    geometric g(0.02); // On average, 2 in 100 products are faulty.
+    double c = 0.95; // 95% confidence.
+    cout << " quantile(g, " << c << ") = " << quantile(g, c) << endl;
+
+    cout << "To be " << c * 100 
+      << "% confident of finding we customer with a fault, need to survey "
+      <<  ceil(quantile(g, c)) << " customers." << endl; // 148
+    c = 0.99; // Very confident.
+    cout << "To be " << c * 100 
+      << "% confident of finding we customer with a fault, need to survey "
+      <<  ceil(quantile(g, c)) << " customers." << endl; // 227
+    c = 0.80; // Only reasonably confident.
+    cout << "To be " << c * 100 
+      << "% confident of finding we customer with a fault, need to survey "
+      <<  ceil(quantile(g, c)) << " customers." << endl; // 79
+
+/*`[h6 Basket Ball Shooters]
+According to Wikipedia, average pro basket ball players get 
+[@http://en.wikipedia.org/wiki/Free_throw free throws]
+in the baskets 70 to 80 % of the time,
+but some get as high as 95%, and others as low as 50%.
+Suppose we want to compare the probabilities
+of failing to get a score only on the first or on the fifth shot?
+To start we will consider the average shooter, say 75%.
+So we construct a geometric distribution
+with success_fraction parameter 75/100 = 0.75.
+*/ 
+    cout.precision(2);
+    geometric gav(0.75); // Shooter averages 7.5 out of 10 in the basket.
+/*`What is probability of getting 1st try in the basket, that is with no failures? */
+    cout << "Probability of score on 1st try = " << pdf(gav, 0) << endl; // 0.75
+/*`This is, of course, the success_fraction probability 75%.
+What is the probability that the shooter only scores on the fifth shot?
+So there are 5-1 = 4 failures before the first success.*/
+    cout << "Probability of score on 5th try = " << pdf(gav, 4) << endl; // 0.0029
+/*`Now compare this with the poor and the best players success fraction.
+We need to constructing new distributions with the different success fractions,
+and then get the corresponding probability density functions values:
+*/
+    geometric gbest(0.95);
+    cout << "Probability of score on 5th try = " << pdf(gbest, 4) << endl; // 5.9e-6
+    geometric gmediocre(0.50);
+    cout << "Probability of score on 5th try = " << pdf(gmediocre, 4) << endl; // 0.031
+/*`So we can see the very much smaller chance (0.000006) of 4 failures by the best shooters,
+compared to the 0.03 of the mediocre.*/
+
+/*`[h6 Estimating failures]
+Of course one man's failure is an other man's success.
+So a fault can be defined as a 'success'.
+
+If a fault occurs once after 100 flights, then one might naively say
+that the risk of fault is obviously 1 in 100 = 1/100, a probability of 0.01.
+
+This is the best estimate we can make, but while it is the truth,
+it is not the whole truth,
+for it hides the big uncertainty when estimating from a single event.
+"One swallow doesn't make a summer."
+To show the magnitude of the uncertainty, the geometric 
+(or the negative binomial) distribution can be used. 
+
+If we chose the popular 95% confidence in the limits, corresponding to an alpha of 0.05,
+because we are calculating a two-sided interval, we must divide alpha by two.
+*/
+    double alpha = 0.05;
+    double k = 100; // So frequency of occurence is 1/100.
+    cout << "Probability is failure is " << 1/k << endl;
+    double t = geometric::find_lower_bound_on_p(k, alpha/2);
+    cout << "geometric::find_lower_bound_on_p(" << int(k) << ", " << alpha/2 << ") = " 
+      << t << endl; // 0.00025
+    t = geometric::find_upper_bound_on_p(k, alpha/2);
+    cout << "geometric::find_upper_bound_on_p(" << int(k) << ", " << alpha/2 << ") = "
+      << t << endl; // 0.037
+/*`So while we estimate the probability is 0.01, it might lie between 0.0003 and 0.04.
+Even if we relax our confidence to alpha = 90%, the bounds only contract to 0.0005 and 0.03.
+And if we require a high confidence, they widen to 0.00005 to 0.05.
+*/
+    alpha = 0.1; // 90% confidence.
+    t = geometric::find_lower_bound_on_p(k, alpha/2);
+    cout << "geometric::find_lower_bound_on_p(" << int(k) << ", " << alpha/2 << ") = "
+      << t << endl; // 0.0005
+    t = geometric::find_upper_bound_on_p(k, alpha/2);
+    cout << "geometric::find_upper_bound_on_p(" << int(k) << ", " << alpha/2 << ") = "
+      << t << endl; // 0.03
+
+    alpha = 0.01; // 99% confidence.
+    t = geometric::find_lower_bound_on_p(k, alpha/2);
+    cout << "geometric::find_lower_bound_on_p(" << int(k) << ", " << alpha/2 << ") = "
+      << t << endl; // 5e-005
+    t = geometric::find_upper_bound_on_p(k, alpha/2);
+    cout << "geometric::find_upper_bound_on_p(" << int(k) << ", " << alpha/2 << ") = "
+        << t << endl; // 0.052
+/*`In real life, there will usually be more than one event (fault or success),
+when the negative binomial, which has the neccessary extra parameter, will be needed.
+*/
+
+/*`As noted above, using a catch block is always a good idea,
+even if you hope not to use it!
+*/
+  }
+  catch(const std::exception& e)
+  { // Since we have set an overflow policy of ignore_error,
+    // an overflow exception should never be thrown.
+     std::cout << "\nMessage from thrown exception was:\n " << e.what() << std::endl;
+/*`
+For example, without a ignore domain error policy, 
+if we asked for ``pdf(g, -1)`` for example, 
+we would get an unhelpful abort, but with a catch:
+[pre
+Message from thrown exception was:
+ Error in function boost::math::pdf(const exponential_distribution<double>&, double):
+ Number of failures argument is -1, but must be >= 0 !
+] [/pre]
+*/
+//] [/ geometric_eg1_2]
+  }
+  return 0;
+}  // int main()
+
+
+/*
+Output is:
+
+  Geometric distribution example
+  
+  success fraction of a six-sided dice is 0.1667
+  0.1667
+  0.1667
+  0.1389
+  pdf(g6, 0) + pdf(g6, 1) = 0.3056
+  cdf(g6, 1) = 0.3056
+  cdf(g6, 12) = 0.9065
+  quantile(g6, 0.99) = 24.26
+  0.1389
+  0.3056
+  quantile(complement(g6, 1 - p))  24.26
+  0.5000346561579232
+  0.5
+   quantile(g, 0.95) = 147.28
+  To be 95% confident of finding we customer with a fault, need to survey 148 customers.
+  To be 99% confident of finding we customer with a fault, need to survey 227 customers.
+  To be 80% confident of finding we customer with a fault, need to survey 79 customers.
+  Probability of score on 1st try = 0.75
+  Probability of score on 5th try = 0.0029
+  Probability of score on 5th try = 5.9e-006
+  Probability of score on 5th try = 0.031
+  Probability is failure is 0.01
+  geometric::find_lower_bound_on_p(100, 0.025) = 0.00025
+  geometric::find_upper_bound_on_p(100, 0.025) = 0.037
+  geometric::find_lower_bound_on_p(100, 0.05) = 0.00051
+  geometric::find_upper_bound_on_p(100, 0.05) = 0.03
+  geometric::find_lower_bound_on_p(100, 0.005) = 5e-005
+  geometric::find_upper_bound_on_p(100, 0.005) = 0.052
+  
+*/
+
+
+
+
+
+
+
+
+
+

Added: trunk/libs/math/example/inverse_chi_squared_example.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/math/example/inverse_chi_squared_example.cpp	2010-11-24 11:34:31 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,215 @@
+// inverse_chi_squared_distribution_example.cpp
+
+// Copyright Paul A. Bristow 2010.
+// Copyright Thomas Mang 2010.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Example 1 of using inverse chi squared distribution
+#include <boost/math/distributions/inverse_chi_squared.hpp>
+using boost::math::inverse_chi_squared_distribution;  // inverse_chi_squared_distribution.
+using boost::math::inverse_chi_squared; //typedef for nverse_chi_squared_distribution double.
+
+#include <iostream>
+using std::cout;    using std::endl;
+#include <iomanip> 
+using std::setprecision;
+using std::setw;
+#include <cmath>
+using std::sqrt;
+
+template <class RealType>
+RealType naive_pdf1(RealType df, RealType x)
+{ // Formula from Wikipedia http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
+  // definition 1 using tgamma for simplicity as a check.
+   using namespace std; // For ADL of std functions.
+   using boost::math::tgamma;
+   RealType df2 = df / 2;
+   RealType result = (pow(2., -df2) * pow(x, (-df2 -1)) * exp(-1/(2 * x) ) )
+      / tgamma(df2);  // 
+   return result;
+}
+
+template <class RealType>
+RealType naive_pdf2(RealType df, RealType x)
+{ // Formula from Wikipedia http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
+  // Definition 2, using tgamma for simplicity as a check.
+  // Not scaled, so assumes scale = 1 and only uses nu aka df.
+   using namespace std; // For ADL of std functions.
+   using boost::math::tgamma;
+   RealType df2 = df / 2;
+   RealType result = (pow(df2, df2) * pow(x, (-df2 -1)) * exp(-df/(2*x) ) )
+     / tgamma(df2);
+   return result;
+}
+
+template <class RealType>
+RealType naive_pdf3(RealType df, RealType scale, RealType x)
+{ // Formula from Wikipedia http://en.wikipedia.org/wiki/Scaled-inverse-chi-square_distribution
+  // *Scaled* version, definition 3, df aka nu, scale aka sigma^2
+  // using tgamma for simplicity as a check.
+   using namespace std; // For ADL of std functions.
+   using boost::math::tgamma;
+   RealType df2 = df / 2;
+   RealType result = (pow(scale * df2, df2) * exp(-df2 * scale/x) ) 
+     / (tgamma(df2) * pow(x, 1+df2));
+   return result;
+}
+
+template <class RealType>
+RealType naive_pdf4(RealType df, RealType scale, RealType x)
+{ // Formula from http://mathworld.wolfram.com/InverseChi-SquaredDistribution.html
+  // Weisstein, Eric W. "Inverse Chi-Squared Distribution." From MathWorld--A Wolfram Web Resource.
+
+  // *Scaled* version, definition 3, df aka nu, scale aka sigma^2
+  // using tgamma for simplicity as a check.
+   using namespace std; // For ADL of std functions.
+   using boost::math::tgamma;
+   RealType nu = df; // Wolfram greek symbols.
+   RealType xi = scale;
+   RealType result = 
+     pow(2, -nu/2) *  exp(- (nu * xi)/(2 * x)) * pow(x, -1-nu/2) * pow((nu * xi), nu/2) 
+     / tgamma(nu/2);
+   return result;
+}
+
+int main()
+{
+
+  cout << "Example (basic) using Inverse chi squared distribution. " << endl;
+
+  cout.precision(std::numeric_limits<double>::max_digits10); // 
+  int i = std::numeric_limits<double>::max_digits10;
+  cout << "Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = " << i << endl; 
+
+  inverse_chi_squared ichsqdef; // All defaults  - not very useful!
+  cout << "default df = " << ichsqdef.degrees_of_freedom()
+    << ", default scale = " <<  ichsqdef.scale() << endl; //  default df = 1, default scale = 0.5
+
+   inverse_chi_squared ichsqdef4(4); // Unscaled version, default scale = 1 / degrees_of_freedom
+   cout << "default df = " << ichsqdef4.degrees_of_freedom()
+    << ", default scale = " <<  ichsqdef4.scale() << endl; //  default df = 4, default scale = 2
+
+   inverse_chi_squared ichsqdef32(3, 2); // Scaled version, both degrees_of_freedom and scale specified.
+   cout << "default df = " << ichsqdef32.degrees_of_freedom()
+    << ", default scale = " <<  ichsqdef32.scale() << endl; // default df = 3, default scale = 2
+
+  {
+    cout.precision(3);
+    double nu = 5.;
+    double scale1 = 1./ nu; // 1st definition sigma^2 = 1/df;
+    double scale2 = 1.; // 2nd definition sigma^2 = 1
+    inverse_chi_squared ichsq(nu, 1/nu); // Not scaled
+    inverse_chi_squared sichsq(nu, 1/nu); // scaled
+
+    cout << "nu = " << ichsq.degrees_of_freedom() << ", scale = " << ichsq.scale() << endl;
+
+    int width = 8;
+
+    cout << " x      pdf(inchsq)    pdf1   pdf2 pdf(scaled)   pdf      pdf    cdf" << endl;
+    for (double x = 0.0; x < 1.; x += 0.1)
+    {
+      cout 
+        << setw(width) << x << ' ' << setw(width) << pdf(ichsq, x) // unscaled
+        << ' ' << setw(width) << naive_pdf1(nu, 1/nu, x) // Wiki def 1 unscaled matches graph 
+        << ' ' << setw(width) << naive_pdf2(nu, scale2, x) // scale = 1 - 2nd definition.
+        << ' ' << setw(width) << naive_pdf3(nu, 1/nu, x) // scaled 
+        << ' ' << setw(width) << naive_pdf4(nu, 1/nu, x) // scaled 
+        << ' ' << setw(width) << pdf(sichsq, x)  // scaled
+        << ' ' << setw(width) << cdf(sichsq, x)  // scaled
+        << ' ' << setw(width) << cdf(ichsq, x) 
+       << endl;
+    }
+  }
+
+  cout.precision(std::numeric_limits<double>::max_digits10);
+
+  inverse_chi_squared ichisq(2., 0.5);
+
+  cout << cdf(ichisq, 1.) << endl;
+
+  return 0;
+}  // int main()
+
+/*
+
+Output is:
+  Example (basic) using Inverse chi squared distribution. 
+  Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = 17
+  default df = 1, default scale = 1
+  default df = 4, default scale = 0.25
+  default df = 3, default scale = 2
+  nu = 5, scale = 0.2
+   x      pdf(inchsq)    pdf1   pdf2 pdf(scaled)   pdf      pdf    cdf
+         0        0    -1.#J    -1.#J    -1.#J    -1.#J        0        0        0
+       0.1     2.83     2.83 3.26e-007     2.83     2.83     2.83   0.0752   0.0752
+       0.2     3.05     3.05  0.00774     3.05     3.05     3.05    0.416    0.416
+       0.3      1.7      1.7    0.121      1.7      1.7      1.7    0.649    0.649
+       0.4    0.941    0.941    0.355    0.941    0.941    0.941    0.776    0.776
+       0.5    0.553    0.553    0.567    0.553    0.553    0.553    0.849    0.849
+       0.6    0.345    0.345    0.689    0.345    0.345    0.345    0.893    0.893
+       0.7    0.227    0.227    0.728    0.227    0.227    0.227    0.921    0.921
+       0.8    0.155    0.155    0.713    0.155    0.155    0.155     0.94     0.94
+       0.9     0.11     0.11    0.668     0.11     0.11     0.11    0.953    0.953
+         1   0.0807   0.0807     0.61   0.0807   0.0807   0.0807    0.963    0.963
+
+*/
+
+
+/*  Parked temporary.
+ double df = 5.;
+  double x = 0.2;
+  // Construction using default RealType double, and default shape and scale..
+  inverse_chi_squared_distribution<> my_inverse_chi_squared(df); // (nu)
+
+  cout << "my_inverse_chi_squared.degrees_of_freedom() = " << my_inverse_chi_squared.degrees_of_freedom() << endl;
+  cout << "x = " << x << ", pdf = " << pdf(my_inverse_chi_squared, x)
+     << ", cdf = " << cdf(my_inverse_chi_squared, x)
+    << endl;
+
+  // Construct using typedef and default shape and scale parameters.
+ 
+  // Example of providing an 'out of domain' or 'bad' parameter,
+  // here a degrees of freedom = 2, for which mean is not defined.
+  // Try block is essential to catch the exception message.
+  // (Uses the default policy which is to throw on all errors).
+  try
+  {
+    inverse_chi_squared ic2(2);
+    cout << "degrees of freedom(ic2) = " << ic2.degrees_of_freedom() << endl;
+
+    cout << "mean(ic2) = " << mean(ic2) << endl;
+  }
+  catch(const std::exception& e)
+  { // Always useful to include try & catch blocks because default policies 
+    // are to throw exceptions on arguments that cause errors like underflow, overflow. 
+    // Lacking try & catch blocks, the program will abort without a message below,
+    // which may give some helpful clues as to the cause of the exception.
+    std::cout <<
+      "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
+  }
+
+  // Example of providing an 'out of domain' or 'bad' parameter,
+  // here a degrees of freedom < 0, for which mean is not defined.
+  // Try block is essential to catch the exception message.
+  // (Uses the default policy which is to throw on all errors).
+  try
+  {
+     inverse_chi_squared icm1(-1);
+
+  }
+  catch(const std::exception& e)
+  { // Always useful to include try & catch blocks because default policies 
+    // are to throw exceptions on arguments that cause errors like underflow, overflow. 
+    // Lacking try & catch blocks, the program will abort without a message below,
+    // which may give some helpful clues as to the cause of the exception.
+    std::cout <<
+      "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
+    //   Message from thrown exception was:
+    // Error in function boost::math::variance(const inverse_gamma_distribution<double>&):
+    // Shape parameter is 1.8999999999999999, but must be > 2!
+  }
+  */
\ No newline at end of file

Added: trunk/libs/math/example/inverse_chi_squared_find_df_example.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/math/example/inverse_chi_squared_find_df_example.cpp	2010-11-24 11:34:31 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,185 @@
+// inverse_chi_squared_distribution_find_df_example.cpp
+
+// Copyright Paul A. Bristow 2010.
+// Copyright Thomas Mang 2010.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+//#define BOOST_MATH_INSTRUMENT
+
+// Example 1 of using inverse chi squared distribution
+#include <boost/math/distributions/inverse_chi_squared.hpp>
+using boost::math::inverse_chi_squared_distribution;  // inverse_chi_squared_distribution.
+using boost::math::inverse_chi_squared; //typedef for nverse_chi_squared_distribution double.
+
+#include <iostream>
+using std::cout;    using std::endl;
+#include <iomanip> 
+using std::setprecision;
+using std::setw;
+#include <cmath>
+using std::sqrt;
+
+int main()
+{
+  cout << "Example using Inverse chi squared distribution to find df. " << endl;
+  try
+  {
+    cout.precision(std::numeric_limits<double>::max_digits10); // 
+    int i = std::numeric_limits<double>::max_digits10;
+    cout << "Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = " << i << endl; 
+
+    cout.precision(3);
+    double nu = 10.;
+    double scale1 = 1./ nu; // 1st definition sigma^2 = 1/df;
+    double scale2 = 1.; // 2nd definition sigma^2 = 1
+    inverse_chi_squared sichsq(nu, 1/nu); // Explicitly scaled to default scale = 1/df.
+    inverse_chi_squared ichsq(nu); // Implicitly scaled to default scale = 1/df.
+    // Try degrees of freedom estimator
+
+    //double df = chi_squared::find_degrees_of_freedom(-diff, alpha[i], alpha[i], variance);
+
+    cout << "ichsq.degrees_of_freedom() = " << ichsq.degrees_of_freedom() << endl;
+
+    double diff = 0.5;  // difference from variance to detect (delta).
+    double variance = 1.; // true variance
+    double alpha = 0.9;
+    double beta = 0.9;
+
+    cout << "diff = " << diff 
+      << ", variance = " << variance << ", ratio = " << diff/variance
+      << ", alpha = " << alpha << ", beta = " << beta << endl;
+    using boost::math::detail::inverse_chi_square_df_estimator;
+    using boost::math::policies::default_policy;
+    inverse_chi_square_df_estimator<> a_df(alpha, beta, variance, diff);
+
+    cout << "df    est" << endl;
+    for (double df = 1; df < 3; df += 0.1)
+    {
+      double est_df = a_df(1);
+      cout << df << "    " << a_df(df) << endl;
+    }
+
+    //template <class F, class T, class Tol, class Policy>std::pair<T, T> 
+    // bracket_and_solve_root(F f, const T& guess, T factor, bool rising, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
+
+
+    //double df = inverse_chi_squared_distribution<>::find_degrees_of_freedom(diff, alpha, beta, variance, 0);
+  
+    double df = inverse_chi_squared::find_degrees_of_freedom(diff, alpha, beta, variance, 100);
+
+    cout << df << endl;
+  }
+  catch(const std::exception& e)
+  { // Always useful to include try & catch blocks because default policies 
+    // are to throw exceptions on arguments that cause errors like underflow, overflow. 
+    // Lacking try & catch blocks, the program will abort without a message below,
+    // which may give some helpful clues as to the cause of the exception.
+    std::cout <<
+      "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
+  }
+  return 0;
+}  // int main()
+
+/*
+
+Output is:
+
+  Example using Inverse chi squared distribution to find df. 
+  Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = 17
+  10
+  
+  Message from thrown exception was:
+     Error in function boost::math::inverse_chi_squared_distribution<double>::inverse_chi_squared_distribution: Degrees of freedom argument is 1.#INF, but must be > 0 !
+diff = 0.5, variance = 1, ratio = 0.5, alpha = 0.1, beta = 0.1
+  df    est
+  1    1
+  ratio+1 = 1.5, quantile(0.1) = 0.00618, cdf = 6.5e-037, result = -0.1
+  1.1    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.00903, cdf = 1.2e-025, result = -0.1
+  1.2    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0125, cdf = 8.25e-019, result = -0.1
+  1.3    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0166, cdf = 2.17e-014, result = -0.1
+  1.4    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0212, cdf = 2.2e-011, result = -0.1
+  1.5    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0265, cdf = 3e-009, result = -0.1
+  1.6    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0323, cdf = 1.11e-007, result = -0.1
+  1.7    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0386, cdf = 1.7e-006, result = -0.1
+  1.8    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0454, cdf = 1.41e-005, result = -0.1
+  1.9    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0527, cdf = 7.55e-005, result = -0.1
+  2    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0604, cdf = 0.000291, result = -0.1
+  2.1    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0685, cdf = 0.00088, result = -0.1
+  2.2    -0.1
+  ratio+1 = 1.5, quantile(0.1) = 0.0771, cdf = 0.0022, result = -0.0999
+  2.3    -0.0999
+  ratio+1 = 1.5, quantile(0.1) = 0.0859, cdf = 0.00475, result = -0.0997
+  2.4    -0.0997
+  ratio+1 = 1.5, quantile(0.1) = 0.0952, cdf = 0.00911, result = -0.0993
+  2.5    -0.0993
+  ratio+1 = 1.5, quantile(0.1) = 0.105, cdf = 0.0159, result = -0.0984
+  2.6    -0.0984
+  ratio+1 = 1.5, quantile(0.1) = 0.115, cdf = 0.0257, result = -0.0967
+  2.7    -0.0967
+  ratio+1 = 1.5, quantile(0.1) = 0.125, cdf = 0.039, result = -0.094
+  2.8    -0.094
+  ratio+1 = 1.5, quantile(0.1) = 0.135, cdf = 0.056, result = -0.0897
+  2.9    -0.0897
+  ratio+1 = 1.5, quantile(0.1) = 20.6, cdf = 1, result = 0.9
+
+    ichsq.degrees_of_freedom() = 10
+  diff = 0.5, variance = 1, ratio = 0.5, alpha = 0.9, beta = 0.9
+  df    est
+  1    1
+  ratio+1 = 1.5, quantile(0.9) = 0.729, cdf = 0.269, result = -0.729
+  1.1    -0.729
+  ratio+1 = 1.5, quantile(0.9) = 0.78, cdf = 0.314, result = -0.693
+  1.2    -0.693
+  ratio+1 = 1.5, quantile(0.9) = 0.83, cdf = 0.36, result = -0.655
+  1.3    -0.655
+  ratio+1 = 1.5, quantile(0.9) = 0.879, cdf = 0.405, result = -0.615
+  1.4    -0.615
+  ratio+1 = 1.5, quantile(0.9) = 0.926, cdf = 0.449, result = -0.575
+  1.5    -0.575
+  ratio+1 = 1.5, quantile(0.9) = 0.973, cdf = 0.492, result = -0.535
+  1.6    -0.535
+  ratio+1 = 1.5, quantile(0.9) = 1.02, cdf = 0.534, result = -0.495
+  1.7    -0.495
+  ratio+1 = 1.5, quantile(0.9) = 1.06, cdf = 0.574, result = -0.455
+  1.8    -0.455
+  ratio+1 = 1.5, quantile(0.9) = 1.11, cdf = 0.612, result = -0.417
+  1.9    -0.417
+  ratio+1 = 1.5, quantile(0.9) = 1.15, cdf = 0.648, result = -0.379
+  2    -0.379
+  ratio+1 = 1.5, quantile(0.9) = 1.19, cdf = 0.681, result = -0.342
+  2.1    -0.342
+  ratio+1 = 1.5, quantile(0.9) = 1.24, cdf = 0.713, result = -0.307
+  2.2    -0.307
+  ratio+1 = 1.5, quantile(0.9) = 1.28, cdf = 0.742, result = -0.274
+  2.3    -0.274
+  ratio+1 = 1.5, quantile(0.9) = 1.32, cdf = 0.769, result = -0.242
+  2.4    -0.242
+  ratio+1 = 1.5, quantile(0.9) = 1.36, cdf = 0.793, result = -0.212
+  2.5    -0.212
+  ratio+1 = 1.5, quantile(0.9) = 1.4, cdf = 0.816, result = -0.184
+  2.6    -0.184
+  ratio+1 = 1.5, quantile(0.9) = 1.44, cdf = 0.836, result = -0.157
+  2.7    -0.157
+  ratio+1 = 1.5, quantile(0.9) = 1.48, cdf = 0.855, result = -0.133
+  2.8    -0.133
+  ratio+1 = 1.5, quantile(0.9) = 1.52, cdf = 0.872, result = -0.11
+  2.9    -0.11
+  ratio+1 = 1.5, quantile(0.9) = 29.6, cdf = 1, result = 0.1
+
+
+  */
\ No newline at end of file

Added: trunk/libs/math/example/inverse_gamma_distribution_example.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/math/example/inverse_gamma_distribution_example.cpp	2010-11-24 11:34:31 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,98 @@
+// inverse_gamma_distribution_example.cpp
+
+// Copyright Paul A. Bristow 2010.
+// Copyright Thomas Mang 2010.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Example 1 of using inverse gamma
+#include <boost/math/distributions/inverse_gamma.hpp>
+using boost::math::inverse_gamma_distribution;  //  inverse_gamma_distribution.
+using boost::math::inverse_gamma;
+
+#include <boost\math\special_functions\gamma.hpp>
+using boost::math::tgamma; // Used for naive pdf as a comparison.
+
+#include <boost/math/distributions/gamma.hpp>
+using boost::math::inverse_gamma_distribution;
+
+#include <iostream>
+using std::cout;    using std::endl;
+#include <iomanip> 
+using std::setprecision;
+#include <cmath>
+using std::sqrt;
+
+int main()
+{
+
+  cout << "Example using Inverse Gamma distribution. " << endl;
+  // TODO - awaiting a real example using Bayesian statistics.
+
+  cout.precision(std::numeric_limits<double>::max_digits10); // 
+
+  int i = std::numeric_limits<double>::max_digits10;
+  cout << "std::numeric_limits<double>::max_digits10 = " << i << endl; 
+
+  double shape = 1.;
+  double scale = 1.;
+  double x = 0.5;
+  // Construction using default RealType double, and default shape and scale..
+  inverse_gamma_distribution<> my_inverse_gamma(shape, scale); // (alpha, beta)
+
+  cout << "my_inverse_gamma.shape() = " << my_inverse_gamma.shape()
+    << ", scale = "<< my_inverse_gamma.scale() << endl;
+  cout << "x = " << x << ", pdf = " << pdf(my_inverse_gamma, x)
+    << ", cdf = " << cdf(my_inverse_gamma, x) << endl;
+
+  // Construct using  typedef and default shape and scale parameters.
+  inverse_gamma my_ig;
+
+  inverse_gamma my_ig23(2, 3);
+  cout << "my_inverse_gamma.shape() = " << my_ig23.shape()
+    << ", scale = "<< my_ig23.scale() << endl;
+  cout << "x = " << x << ", pdf = " << pdf(my_ig23, x)
+    << ", cdf = " << cdf(my_ig23, x) << endl;
+
+  // Example of providing an 'out of domain' or 'bad' parameter,
+  // here a shape < 1, for which mean is not defined.
+  // Try block is essential to catch the exception message.
+  // (Uses the default policy which is to throw on all errors).
+  try
+  {
+    inverse_gamma if051(0.5, 1);
+    //inverse_gamma if051(0.5, 1);
+    cout << "mean(if051) = " << mean(if051) << endl;
+  }
+  catch(const std::exception& e)
+  { // Always useful to include try & catch blocks because default policies 
+    // are to throw exceptions on arguments that cause errors like underflow, overflow. 
+    // Lacking try & catch blocks, the program will abort without a message below,
+    // which may give some helpful clues as to the cause of the exception.
+    std::cout <<
+      "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
+  }
+
+  return 0;
+}  // int main()
+
+/*
+
+Output is:
+  Example using Inverse Gamma distribution. 
+  std::numeric_limits<double>::max_digits10 = 17
+  my_inverse_gamma.shape() = 1, scale = 1
+  x = 0.5, pdf = 0.54134113294645081, cdf = 0.1353352832366127
+  my_inverse_gamma.shape() = 2, scale = 3
+  x = 0.5, pdf = 0.17847015671997774, cdf = 0.017351265236664509
+  
+  Message from thrown exception was:
+     Error in function boost::math::mean(const inverse_gamma_distribution<double>&): Shape parameter is 0.5, but for a defined mean it must be > 1
+
+
+*/
+
+

Added: trunk/libs/math/example/inverse_gamma_example.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/math/example/inverse_gamma_example.cpp	2010-11-24 11:34:31 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,49 @@
+// inverse_gamma_example.cpp
+
+// Copyright Paul A. Bristow 2010.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Example 1 of using inverse gamma functions.
+
+#include <boost/math/special_functions/gamma.hpp>
+
+using boost::math::gamma_p_inv; // Compute x given a
+//using boost::math::gamma_q_inv;
+//using boost::math::gamma_p_inva; // Compute a given x
+//using boost::math::gamma_q_inva;
+
+#include <iostream>
+   using std::cout;    using std::endl;
+#include <iomanip>
+   using std::setprecision;
+#include <cmath>
+   using std::sqrt;
+#include <limits>
+
+int main()
+{
+  cout << "Example 1 using Inverse Gamma function. " << endl;
+
+  cout.precision(std::numeric_limits<double>::max_digits10);
+  
+  double x = 1.;
+  double a = 10;
+
+  double r = boost::math::gamma_q_inv(a ,x);
+
+  cout << " x = " << x << ", = gamma_q_inv(a,x)" << r << endl; // 
+  
+  return 0;
+}  // int main()
+
+/*
+
+Output is:
+
+*/
+
+