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