]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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> | |
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> | |
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> | |
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> | |
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 | ||
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 |