
BoostCommit : 
Subject: [Boostcommit] svn:boost r73047  trunk/libs/math/example
From: pbristow_at_[hidden]
Date: 20110713 10:19:35
Author: pbristow
Date: 20110713 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 20110713 10:19:33 EDT (Wed, 13 Jul 2011)
@@ 39,14 +39,14 @@
//[inverse_chi_squared_bayes_eg_1
/*`
The scaledinversedchisquared 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 scaledinversechisquared 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 scaledinversechisquared distributed.
@@ 95,23 +95,24 @@
Using the survey results, we try to find the parameter set
of a scaledinversechisquared 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 scaledinversechisquared distribution
as guideline, and using some trialanderror 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
chisquared 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 scaledinversechisquared 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 trialanderror and calls to the global quantile function, we also find that a
+value of 20 for the degreesoffreedom ([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 onehalf (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
+onehalf (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 ['mustnot] 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 ['mustnot]) 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 scaledinversechisquare distributed,
then the posterior distribution is also scaledinversechisquareddistributed
(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 chisquare 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.
Machinespecific 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,
+Machinespecific 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 scaledinversechisquared distribution with df = 120 and scale = 49.54.
+is scaledinversechisquared distribution with degreesoffreedom ([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 degreesoffreedom: " << priorDF << endl;
+ cout << "posterior degreesoffreedom: " << 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 degreesoffreedom: 20
+ posterior degreesoffreedom: 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 onemillionth of a permille).
+and find a distribution clearly shifted to greater values (see figure).
On the other hand, with an almost nearhalf 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 onemillionth of a permille).
+On the other hand, with an almost nearhalf 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 degreesoffreedom: 20
+ posterior degreesoffreedom: 120
+ prior scale: 25
+ posterior scale: 49.5
Posterior distribution:
2.5% quantile: 39.1
BoostCommit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk