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

Subject: [Boost-commit] svn:boost r73047 - trunk/libs/math/example
From: pbristow_at_[hidden]
Date: 2011-07-13 10:19:35

Author: pbristow
Date: 2011-07-13 10:19:33 EDT (Wed, 13 Jul 2011)
New Revision: 73047
URL: http://svn.boost.org/trac/boost/changeset/73047

Log:
New Bayes example
Text files modified: 
   trunk/libs/math/example/inverse_chi_squared_bayes_eg.cpp |   145 +++++++++++++++++++++------------------ 
   1 files changed, 79 insertions(+), 66 deletions(-)

Modified: trunk/libs/math/example/inverse_chi_squared_bayes_eg.cpp
==============================================================================
--- trunk/libs/math/example/inverse_chi_squared_bayes_eg.cpp	(original)
+++ trunk/libs/math/example/inverse_chi_squared_bayes_eg.cpp	2011-07-13 10:19:33 EDT (Wed, 13 Jul 2011)
@@ -39,14 +39,14 @@
 //[inverse_chi_squared_bayes_eg_1
 /*`
 The scaled-inversed-chi-squared distribution is the conjugate prior distribution
-for the variance ([sigma] [super 2]) parameter of a normal distribution
+for the variance ([sigma][super 2]) parameter of a normal distribution
 with known expectation ([mu]).
 As such it has widespread application in Bayesian statistics:
 
 In [@http://en.wikipedia.org/wiki/Bayesian_inference Bayesian inference],
 the strength of belief into certain parameter values is
-itself described through a distribution, and parameters
-hence become themselves random variables.
+itself described through a distribution. Parameters
+hence become themselves modelled and interpreted as random variables.
 
 In this worked example, we perform such a Bayesian analysis by using
 the scaled-inverse-chi-squared distribution as prior and posterior distribution
@@ -67,19 +67,19 @@
 Ideally each ball should have a diameter of precisely 3000 [mu]m (3 mm).
 Assume that machines generally produce balls of that size on average (mean),
 but individual balls can vary slightly in either direction
-following a normal distribution. Depending on various production conditions
+following (approximately) a normal distribution. Depending on various production conditions
 (e.g. raw material used for balls, workplace temperature and humidity, maintenance frequency and quality)
 some machines produce balls tighter distributed around the target of 3000 [mu]m,
 while others produce balls with a wider distribution. 
 Therefore the variance parameter of the normal distribution of the ball sizes varies
 from machine to machine. An extensive survey by the precision machinery manufacturer, however,
 has shown that most machines operate with a variance between 15 and 50,
-and near 25 [mu]m[super 2] on average (or to a standard deviation of 5 [mu]m).
+and near 25 [mu]m[super 2] on average.
 
 Using this information, we want to model the variance of the machines.
 The variance is strictly positive, and therefore we look for a statistical distribution
 with support in the positive domain of the real numbers.
-Because the expectation of the normal distribution of the balls is known (3000 [mu]m),
+Given the expectation of the normal distribution of the balls is known (3000 [mu]m),
 for reasons of conjugacy, it is customary practice in Bayesian statistics
 to model the variance to be scaled-inverse-chi-squared distributed. 
 
@@ -95,23 +95,24 @@
 
 Using the survey results, we try to find the parameter set
 of a scaled-inverse-chi-squared distribution 
-so that the properties of this distribution match the results.
+so that the properties of this distribution match the results. 
 Using the mathematical properties of the scaled-inverse-chi-squared distribution 
-as guideline, and using some trial-and-error and calls to the global quantile function,
-we find the set of degrees of freedom df = 20 and scale = 25 [mu]m to be adequate.
-These values df = 20 and scale = 25 are chosen to get a scaled inverse 
-chi-squared distribution having the properties of most mass between 
-15 and 50, and the mode near 25.
+as guideline, we see that that both the mean and mode of the scaled-inverse-chi-squared distribution
+are approximately given by the scale parameter (s) of the distribution. As the survey machines operated at a
+variance of 25 [mu]m[super 2] on average, we hence set the scale parameter (s[sub prior]) of our prior distribution 
+equal to this value. Using some trial-and-error and calls to the global quantile function, we also find that a
+value of 20 for the degrees-of-freedom ([nu][sub prior]) parameter is adequate so that
+most of the prior distribution mass is located between 15 and 50 (see figure below).
 
-We first construct our prior using these values, and then list out a few quantiles.
+We first construct our prior distribution using these values, and then list out a few quantiles:
 
 */
   double priorDF = 20.0;
-  double priorScale =  25.0; 
+  double priorScale = 25.0; 
 
   inverse_chi_squared prior(priorDF, priorScale);
-  // Using the inverse_gamma distribution instead, we coudl equivalently write
-  // inverse_gamma prior(priorDF /2, priorScale * priorDF /2);
+  // Using an inverse_gamma distribution instead, we could equivalently write
+  // inverse_gamma prior(priorDF / 2.0, priorScale * priorDF / 2.0);
   
   cout << "Prior distribution:" << endl << endl;
   cout << "  2.5% quantile: " << quantile(prior, 0.025) << endl;
@@ -135,7 +136,7 @@
 //[inverse_chi_squared_bayes_eg_2
 /*`
 Based on this distribution, we can now calculate the probability of having a machine
-with precision <= 15 or > 50.
+working with an unusual work precision (variance) at <= 15 or > 50.
 For this task, we use calls to the `boost::math::` functions `cdf` and `complement`,
 respectively, and find a probability of about 0.031 (3.1%) for each case.
 */
@@ -160,86 +161,87 @@
 //[inverse_chi_squared_bayes_eg_3
 /*`Therefore, only 3.1% of all precision machines produce balls with a variance of 15 or less
 (particularly precise machines),
-but also only 3.1% of all machines produce balls
-with a variance of as high as 50 or more (particularly imprecise machines).
-Only with a probability of just over one-half (1 - 0.45793 = 54.2%)
-is the variance actually less than 25.
+but also only 3.2% of all machines produce balls
+with a variance of as high as 50 or more (particularly imprecise machines). Moreover, slightly more than
+one-half (1 - 0.458 = 54.2%) of the machines work at a variance greater than 25. 
 
 Notice here the distinction between a
 [@http://en.wikipedia.org/wiki/Bayesian_inference Bayesian] analysis and a
 [@http://en.wikipedia.org/wiki/Frequentist_inference frequentist] analysis:
 because we model the variance as random variable itself,
-we can calculate and straightforwardly interpret probabilities for the parameter directly,
-which is generally a strict ['must-not] in the frequentist world.
+we can calculate and straightforwardly interpret probabilities for given parameter values directly,
+while such an approach is not possible (and interpretationally a strict ['must-not]) in the frequentist
+world.
 
 [h5 Step 2: Investigate a single machine]
 
 In the second step, we investigate a single machine,
 which is suspected to suffer from a major fault
 as the produced balls show fairly high size variability.
-Based on the prior distribution of generic machinery performed (derived above)
-and data on produced balls, we calculate the posterior distribution for that machine
-and use its properties for guidance regarding continued machine operation or suspension.
+Based on the prior distribution of generic machinery performance (derived above)
+and data on balls produced by the suspect machine, we calculate the posterior distribution for that 
+machine and use its properties for guidance regarding continued machine operation or suspension.
 
 It can be shown that if the prior distribution
 was chosen to be scaled-inverse-chi-square distributed,
 then the posterior distribution is also scaled-inverse-chi-squared-distributed 
 (prior and posterior distributions are hence conjugate).
-For more details regarding conjugates and formula to derive the parameters set
+For more details regarding conjugacy and formula to derive the parameters set
 for the posterior distribution see
 [@http://en.wikipedia.org/wiki/Conjugate_prior Conjugate prior].
 
-From table of conjugate distributions, the Posterior hyperparameters are given,
-and for the Scaled inverse chi-square distribution
-(and normal likelihood and known mean, and model parameter variance)
-is given by the two expressions given for posterior degrees of freedom and scale:
 
-__spaces [nu] = [nu] + n
+Given the prior distribution parameters and sample data (of size n), the posterior distribution parameters 
+are given by the two expressions:
+
+__spaces [nu][sub posterior] = [nu][sub prior] + n
 
 which gives the posteriorDF below, and
 
-__spaces [sigma][super 2] = [nu][sigma][super 2] + [Sigma][super n][sub i=1](x[sub i] - [mu])[super 2]/ ([nu] + n)
+__spaces s[sub posterior] = ([nu][sub prior]s[sub prior] + [Sigma][super n][sub i=1](x[sub i] - [mu])[super 2]) / ([nu][sub prior] + n)
 
 which after some rearrangement gives the formula for the posteriorScale below.
 
-Machine-specific data consists of 100 balls which were accurately measured
-and show a mean of 25 [mu]m and a variance of 55.
-From this data, and the prior parameterization,
+Machine-specific data consist of 100 balls which were accurately measured
+and show the expected mean of 3000 [mu]m and a sample variance of 55 (calculated for a sample mean defined to be 3000 exactly).
+From these data, the prior parameterization, and noting that the term 
+[Sigma][super n][sub i=1](x[sub i] - [mu])[super 2] equals the sample variance multiplied by n - 1,
 it follows that the posterior distribution of the variance parameter
-is scaled-inverse-chi-squared distribution with df = 120 and scale = 49.54.
+is scaled-inverse-chi-squared distribution with degrees-of-freedom ([nu][sub posterior]) = 120 and 
+scale (s[sub posterior]) = 49.54.
 */
 
   int ballsSampleSize = 100;
-  cout <<"balls sample Size " << ballsSampleSize << endl;
+  cout <<"balls sample size: " << ballsSampleSize << endl;
   double ballsSampleVariance = 55.0;
-  cout <<"balls sample variance " << ballsSampleVariance << endl;
+  cout <<"balls sample variance: " << ballsSampleVariance << endl;
 
   double posteriorDF = priorDF + ballsSampleSize;
-  cout << "prior degrees of freedom " << priorDF << endl;
-  cout << "Posterior degrees of freedom " << posteriorDF << endl;
+  cout << "prior degrees-of-freedom: " << priorDF << endl;
+  cout << "posterior degrees-of-freedom: " << posteriorDF << endl;
   
   double posteriorScale = 
     (priorDF * priorScale + (ballsSampleVariance * (ballsSampleSize - 1))) / posteriorDF;
-  cout << "Prior scale " << priorScale  << endl;
-  cout << "Posterior scale " << posteriorScale << endl;
+  cout << "prior scale: " << priorScale  << endl;
+  cout << "posterior scale: " << posteriorScale << endl;
 
 /*`An interesting feature here is that one needs only to know a summary statistics of the sample
-to specify the posterior parameters: the 100 individual measurements are irrelevant,
-just knowledge of the variance and number of measurements is sufficient.
+to parameterize the posterior distribution: the 100 individual ball measurements are irrelevant,
+just knowledge of the sample variance and number of measurements is sufficient.
 */
 
 //] [/inverse_chi_squared_bayes_eg_3]
 
 //[inverse_chi_squared_bayes_eg_output_3
-/*`that produces this output:
+/*`That produces this output:
 
 
-  balls sample Size 100
-  balls sample variance 55
-  prior degrees of freedom 20
-  Posterior degrees of freedom 120
-  Prior scale 25
-  Posterior scale 49.5
+  balls sample size: 100
+  balls sample variance: 55
+  prior degrees-of-freedom: 20
+  posterior degrees-of-freedom: 120
+  prior scale: 25
+  posterior scale: 49.5
    
   */
 //] [/inverse_chi_squared_bayes_eg_output_3]
@@ -247,18 +249,13 @@
 //[inverse_chi_squared_bayes_eg_4
 /*`To compare the generic machinery performance with our suspect machine,
 we calculate again the same quantiles and probabilities as above,
-and find a distribution clearly shifted to greater values.
-Indeed, the probability that the machine works at a low variance (<= 15) is almost zero,
-and even the probability of working at average or better performance is negligibly small
-(less than one-millionth of a permille).
+and find a distribution clearly shifted to greater values (see figure).
 
-On the other hand, with an almost near-half probability (49%), the variance is actually
-in the extreme high variance range of > 50.
+[graph prior_posterior_plot]
 
-Based on this information the operation of the machine is taken out of use and serviced.
 */
 
-  inverse_chi_squared posterior(posteriorDF, posteriorScale);
+ inverse_chi_squared posterior(posteriorDF, posteriorScale);
 
   cout << "Posterior distribution:" << endl << endl;
   cout << "  2.5% quantile: " << boost::math::quantile(posterior, 0.025) << endl;
@@ -269,6 +266,7 @@
   cout << "  probability variance <= 25: " << boost::math::cdf(posterior, 25.0) << endl;
   cout << "  probability variance > 50: " 
     << boost::math::cdf(boost::math::complement(posterior, 50.0)) << endl;
+
 //] [/inverse_chi_squared_bayes_eg_4]
 
 //[inverse_chi_squared_bayes_eg_output_4
@@ -287,6 +285,21 @@
   */
 //] [/inverse_chi_squared_bayes_eg_output_4]
 
+//[inverse_chi_squared_bayes_eg_5
+/*`Indeed, the probability that the machine works at a low variance (<= 15) is almost zero,
+and even the probability of working at average or better performance is negligibly small
+(less than one-millionth of a permille). 
+On the other hand, with an almost near-half probability (49%), the machine operates in the
+extreme high variance range of > 50 characteristic for poorly performing machines.
+
+Based on this information the operation of the machine is taken out of use and serviced.
+
+In summary, the Bayesian analysis allowed us to make exact probabilistic statements about a
+parameter of interest, and hence provided us results with straightforward interpretation. 
+
+*/
+//] [/inverse_chi_squared_bayes_eg_5]
+
 } // int main()
 
 //[inverse_chi_squared_bayes_eg_output
@@ -304,12 +317,12 @@
     probability variance <= 25: 0.458
     probability variance > 50: 0.0318
   
-  balls sample Size 100
-  balls sample variance 55
-  prior degrees of freedom 20
-  Posterior degrees of freedom 120
-  Prior scale 25
-  Posterior scale 49.5
+  balls sample size: 100
+  balls sample variance: 55
+  prior degrees-of-freedom: 20
+  posterior degrees-of-freedom: 120
+  prior scale: 25
+  posterior scale: 49.5
   Posterior distribution:
   
     2.5% quantile: 39.1