| // inverse_chi_squared_bayes_eg.cpp |
| |
| // Copyright Thomas Mang 2011. |
| // Copyright Paul A. Bristow 2011. |
| |
| // 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! |
| |
| #include <iostream> |
| // using std::cout; using std::endl; |
| |
| //#define define possible error-handling macros here? |
| |
| #include "boost/math/distributions.hpp" |
| // using ::boost::math::inverse_chi_squared; |
| |
| int main() |
| { |
| using std::cout; using std::endl; |
| |
| using ::boost::math::inverse_chi_squared; |
| using ::boost::math::inverse_gamma; |
| using ::boost::math::quantile; |
| using ::boost::math::cdf; |
| |
| cout << "Inverse_chi_squared_distribution Bayes example: " << endl <<endl; |
| |
| cout.precision(3); |
| // Examples of using the inverse_chi_squared distribution. |
| |
| //[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 |
| 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. 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 |
| for the variance parameter of a normal distribution. |
| |
| For more general information on Bayesian type of analyses, |
| see: |
| |
| * Andrew Gelman, John B. Carlin, Hal E. Stern, Donald B. Rubin, Bayesian Data Analysis, |
| 2003, ISBN 978-1439840955. |
| |
| * Jim Albert, Bayesian Compution with R, Springer, 2009, ISBN 978-0387922973. |
| |
| (As the scaled-inversed-chi-squared is another parameterization of the inverse-gamma distribution, |
| this example could also have used the inverse-gamma distribution). |
| |
| Consider precision machines which produce balls for a high-quality ball bearing. |
| 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 (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. |
| |
| 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. |
| 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. |
| |
| In a first step, we will try to use the survey information to model |
| the general knowledge about the variance parameter of machines measured by the manufacturer. |
| This will provide us with a generic prior distribution that is applicable |
| if nothing more specific is known about a particular machine. |
| |
| In a second step, we will then combine the prior-distribution information in a Bayesian analysis |
| with data on a specific single machine to derive a posterior distribution for that machine. |
| |
| [h5 Step one: Using the survey information.] |
| |
| 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. |
| Using the mathematical properties of the scaled-inverse-chi-squared distribution |
| 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 distribution using these values, and then list out a few quantiles: |
| |
| */ |
| double priorDF = 20.0; |
| double priorScale = 25.0; |
| |
| inverse_chi_squared prior(priorDF, priorScale); |
| // 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; |
| cout << " 50% quantile: " << quantile(prior, 0.5) << endl; |
| cout << " 97.5% quantile: " << quantile(prior, 0.975) << endl << endl; |
| |
| //] [/inverse_chi_squared_bayes_eg_1] |
| |
| //[inverse_chi_squared_bayes_eg_output_1 |
| /*`This produces this output: |
| |
| Prior distribution: |
| |
| 2.5% quantile: 14.6 |
| 50% quantile: 25.9 |
| 97.5% quantile: 52.1 |
| |
| */ |
| //] [/inverse_chi_squared_bayes_eg_output_1] |
| |
| //[inverse_chi_squared_bayes_eg_2 |
| /*` |
| Based on this distribution, we can now calculate the probability of having a machine |
| 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. |
| */ |
| |
| cout << " probability variance <= 15: " << boost::math::cdf(prior, 15.0) << endl; |
| cout << " probability variance <= 25: " << boost::math::cdf(prior, 25.0) << endl; |
| cout << " probability variance > 50: " |
| << boost::math::cdf(boost::math::complement(prior, 50.0)) |
| << endl << endl; |
| //] [/inverse_chi_squared_bayes_eg_2] |
| |
| //[inverse_chi_squared_bayes_eg_output_2 |
| /*`This produces this output: |
| |
| probability variance <= 15: 0.031 |
| probability variance <= 25: 0.458 |
| probability variance > 50: 0.0318 |
| |
| */ |
| //] [/inverse_chi_squared_bayes_eg_output_2] |
| |
| //[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.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 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 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 conjugacy and formula to derive the parameters set |
| for the posterior distribution see |
| [@http://en.wikipedia.org/wiki/Conjugate_prior Conjugate prior]. |
| |
| |
| 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 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 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 degrees-of-freedom ([nu][sub posterior]) = 120 and |
| scale (s[sub posterior]) = 49.54. |
| */ |
| |
| int ballsSampleSize = 100; |
| cout <<"balls sample size: " << ballsSampleSize << endl; |
| double ballsSampleVariance = 55.0; |
| cout <<"balls sample variance: " << ballsSampleVariance << endl; |
| |
| double posteriorDF = priorDF + ballsSampleSize; |
| 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; |
| |
| /*`An interesting feature here is that one needs only to know a summary statistics of the sample |
| 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: |
| |
| |
| 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] |
| |
| //[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 (see figure). |
| |
| [graph prior_posterior_plot] |
| |
| */ |
| |
| inverse_chi_squared posterior(posteriorDF, posteriorScale); |
| |
| cout << "Posterior distribution:" << endl << endl; |
| cout << " 2.5% quantile: " << boost::math::quantile(posterior, 0.025) << endl; |
| cout << " 50% quantile: " << boost::math::quantile(posterior, 0.5) << endl; |
| cout << " 97.5% quantile: " << boost::math::quantile(posterior, 0.975) << endl << endl; |
| |
| cout << " probability variance <= 15: " << boost::math::cdf(posterior, 15.0) << endl; |
| 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 |
| /*`This produces this output: |
| |
| Posterior distribution: |
| |
| 2.5% quantile: 39.1 |
| 50% quantile: 49.8 |
| 97.5% quantile: 64.9 |
| |
| probability variance <= 15: 2.97e-031 |
| probability variance <= 25: 8.85e-010 |
| probability variance > 50: 0.489 |
| |
| */ |
| //] [/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 |
| /*` |
| [pre |
| Inverse_chi_squared_distribution Bayes example: |
| |
| Prior distribution: |
| |
| 2.5% quantile: 14.6 |
| 50% quantile: 25.9 |
| 97.5% quantile: 52.1 |
| |
| probability variance <= 15: 0.031 |
| 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 |
| Posterior distribution: |
| |
| 2.5% quantile: 39.1 |
| 50% quantile: 49.8 |
| 97.5% quantile: 64.9 |
| |
| probability variance <= 15: 2.97e-031 |
| probability variance <= 25: 8.85e-010 |
| probability variance > 50: 0.489 |
| |
| ] [/pre] |
| */ |
| //] [/inverse_chi_squared_bayes_eg_output] |