1 // inverse_chi_squared_distribution_example.cpp
2
3 // Copyright Paul A. Bristow 2010.
4 // Copyright Thomas Mang 2010.
5
6 // Use, modification and distribution are subject to the
7 // Boost Software License, Version 1.0.
8 // (See accompanying file LICENSE_1_0.txt
9 // or copy at http://www.boost.org/LICENSE_1_0.txt)
10
11 // Example 1 of using inverse chi squared distribution
12 #include <boost/math/distributions/inverse_chi_squared.hpp>
13 using boost::math::inverse_chi_squared_distribution; // inverse_chi_squared_distribution.
14 using boost::math::inverse_chi_squared; //typedef for nverse_chi_squared_distribution double.
15
16 #include <iostream>
17 using std::cout; using std::endl;
18 #include <iomanip>
19 using std::setprecision;
20 using std::setw;
21 #include <cmath>
22 using std::sqrt;
23
24 template <class RealType>
naive_pdf1(RealType df,RealType x)25 RealType naive_pdf1(RealType df, RealType x)
26 { // Formula from Wikipedia http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
27 // definition 1 using tgamma for simplicity as a check.
28 using namespace std; // For ADL of std functions.
29 using boost::math::tgamma;
30 RealType df2 = df / 2;
31 RealType result = (pow(2., -df2) * pow(x, (-df2 -1)) * exp(-1/(2 * x) ) )
32 / tgamma(df2); //
33 return result;
34 }
35
36 template <class RealType>
naive_pdf2(RealType df,RealType x)37 RealType naive_pdf2(RealType df, RealType x)
38 { // Formula from Wikipedia http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
39 // Definition 2, using tgamma for simplicity as a check.
40 // Not scaled, so assumes scale = 1 and only uses nu aka df.
41 using namespace std; // For ADL of std functions.
42 using boost::math::tgamma;
43 RealType df2 = df / 2;
44 RealType result = (pow(df2, df2) * pow(x, (-df2 -1)) * exp(-df/(2*x) ) )
45 / tgamma(df2);
46 return result;
47 }
48
49 template <class RealType>
naive_pdf3(RealType df,RealType scale,RealType x)50 RealType naive_pdf3(RealType df, RealType scale, RealType x)
51 { // Formula from Wikipedia http://en.wikipedia.org/wiki/Scaled-inverse-chi-square_distribution
52 // *Scaled* version, definition 3, df aka nu, scale aka sigma^2
53 // using tgamma for simplicity as a check.
54 using namespace std; // For ADL of std functions.
55 using boost::math::tgamma;
56 RealType df2 = df / 2;
57 RealType result = (pow(scale * df2, df2) * exp(-df2 * scale/x) )
58 / (tgamma(df2) * pow(x, 1+df2));
59 return result;
60 }
61
62 template <class RealType>
naive_pdf4(RealType df,RealType scale,RealType x)63 RealType naive_pdf4(RealType df, RealType scale, RealType x)
64 { // Formula from http://mathworld.wolfram.com/InverseChi-SquaredDistribution.html
65 // Weisstein, Eric W. "Inverse Chi-Squared Distribution." From MathWorld--A Wolfram Web Resource.
66 // *Scaled* version, definition 3, df aka nu, scale aka sigma^2
67 // using tgamma for simplicity as a check.
68 using namespace std; // For ADL of std functions.
69 using boost::math::tgamma;
70 RealType nu = df; // Wolfram uses greek symbols nu,
71 RealType xi = scale; // and xi.
72 RealType result =
73 pow(2, -nu/2) * exp(- (nu * xi)/(2 * x)) * pow(x, -1-nu/2) * pow((nu * xi), nu/2)
74 / tgamma(nu/2);
75 return result;
76 }
77
main()78 int main()
79 {
80
81 cout << "Example (basic) using Inverse chi squared distribution. " << endl;
82
83 // TODO find a more practical/useful example. Suggestions welcome?
84
85 #ifdef BOOST_NO_CXX11_NUMERIC_LIMITS
86 int max_digits10 = 2 + (boost::math::policies::digits<double, boost::math::policies::policy<> >() * 30103UL) / 100000UL;
87 cout << "BOOST_NO_CXX11_NUMERIC_LIMITS is defined" << endl;
88 #else
89 int max_digits10 = std::numeric_limits<double>::max_digits10;
90 #endif
91 cout << "Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = "
92 << max_digits10 << endl;
93 cout.precision(max_digits10); //
94
95 inverse_chi_squared ichsqdef; // All defaults - not very useful!
96 cout << "default df = " << ichsqdef.degrees_of_freedom()
97 << ", default scale = " << ichsqdef.scale() << endl; // default df = 1, default scale = 0.5
98
99 inverse_chi_squared ichsqdef4(4); // Unscaled version, default scale = 1 / degrees_of_freedom
100 cout << "default df = " << ichsqdef4.degrees_of_freedom()
101 << ", default scale = " << ichsqdef4.scale() << endl; // default df = 4, default scale = 2
102
103 inverse_chi_squared ichsqdef32(3, 2); // Scaled version, both degrees_of_freedom and scale specified.
104 cout << "default df = " << ichsqdef32.degrees_of_freedom()
105 << ", default scale = " << ichsqdef32.scale() << endl; // default df = 3, default scale = 2
106
107 {
108 cout.precision(3);
109 double nu = 5.;
110 //double scale1 = 1./ nu; // 1st definition sigma^2 = 1/df;
111 //double scale2 = 1.; // 2nd definition sigma^2 = 1
112 inverse_chi_squared ichsq(nu, 1/nu); // Not scaled
113 inverse_chi_squared sichsq(nu, 1/nu); // scaled
114
115 cout << "nu = " << ichsq.degrees_of_freedom() << ", scale = " << ichsq.scale() << endl;
116
117 int width = 8;
118
119 cout << " x pdf pdf1 pdf2 pdf(scaled) pdf pdf cdf cdf" << endl;
120 for (double x = 0.0; x < 1.; x += 0.1)
121 {
122 cout
123 << setw(width) << x
124 << ' ' << setw(width) << pdf(ichsq, x) // unscaled
125 << ' ' << setw(width) << naive_pdf1(nu, x) // Wiki def 1 unscaled matches graph
126 << ' ' << setw(width) << naive_pdf2(nu, x) // scale = 1 - 2nd definition.
127 << ' ' << setw(width) << naive_pdf3(nu, 1/nu, x) // scaled
128 << ' ' << setw(width) << naive_pdf4(nu, 1/nu, x) // scaled
129 << ' ' << setw(width) << pdf(sichsq, x) // scaled
130 << ' ' << setw(width) << cdf(sichsq, x) // scaled
131 << ' ' << setw(width) << cdf(ichsq, x) // unscaled
132 << endl;
133 }
134 }
135
136 cout.precision(max_digits10);
137
138 inverse_chi_squared ichisq(2., 0.5);
139 cout << "pdf(ichisq, 1.) = " << pdf(ichisq, 1.) << endl;
140 cout << "cdf(ichisq, 1.) = " << cdf(ichisq, 1.) << endl;
141
142 return 0;
143 } // int main()
144
145 /*
146
147 Output is:
148 Example (basic) using Inverse chi squared distribution.
149 Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = 17
150 default df = 1, default scale = 1
151 default df = 4, default scale = 0.25
152 default df = 3, default scale = 2
153 nu = 5, scale = 0.2
154 x pdf pdf1 pdf2 pdf(scaled) pdf pdf cdf cdf
155 0 0 -1.#J -1.#J -1.#J -1.#J 0 0 0
156 0.1 2.83 2.83 3.26e-007 2.83 2.83 2.83 0.0752 0.0752
157 0.2 3.05 3.05 0.00774 3.05 3.05 3.05 0.416 0.416
158 0.3 1.7 1.7 0.121 1.7 1.7 1.7 0.649 0.649
159 0.4 0.941 0.941 0.355 0.941 0.941 0.941 0.776 0.776
160 0.5 0.553 0.553 0.567 0.553 0.553 0.553 0.849 0.849
161 0.6 0.345 0.345 0.689 0.345 0.345 0.345 0.893 0.893
162 0.7 0.227 0.227 0.728 0.227 0.227 0.227 0.921 0.921
163 0.8 0.155 0.155 0.713 0.155 0.155 0.155 0.94 0.94
164 0.9 0.11 0.11 0.668 0.11 0.11 0.11 0.953 0.953
165 1 0.0807 0.0807 0.61 0.0807 0.0807 0.0807 0.963 0.963
166 pdf(ichisq, 1.) = 0.30326532985631671
167 cdf(ichisq, 1.) = 0.60653065971263365
168
169
170 */
171
172