-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathchisquared.h
More file actions
192 lines (149 loc) · 4.98 KB
/
chisquared.h
File metadata and controls
192 lines (149 loc) · 4.98 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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
#ifndef CHISQUARED_H
#define CHISQUARED_H
#include "statistics/gamma.h"
#include "statistics/normal.h"
#include "statistics/prob.h"
#include <cassert>
#include <cmath>
namespace Stats {
//! Chi-squared distribution
class ChiSquared
{
public:
explicit ChiSquared( int nu ); //!< Value Constructor
ChiSquared (const ChiSquared & other) = delete; //!< Copy constructor
ChiSquared (const ChiSquared && other) = delete; //!< Move constructor
ChiSquared & operator = (const ChiSquared & other) = delete; //!< Assignment operator
ChiSquared & operator = (const ChiSquared && other) = delete; //!< Move assignment operator
~ChiSquared() = default;
[[nodiscard]] double pdf( double x) const; //!< Probability density function
[[nodiscard]] Prob cdf( double x) const; //!< Cumulative distribution function
[[nodiscard]] double icdf( Prob P) const; //!< Inverse cumulative distribution function
[[nodiscard]] double mean() const { return nu; } //!< Mean
[[nodiscard]] double var() const { return 2*nu; } //!< Variance
[[nodiscard]] double mode() const { return std::fmax( nu-2, 0.); } //!< mode
[[nodiscard]] double rnd() const; //!< Random number
[[nodiscard]] int dof() const { return nu; } //!< Get degrees of freedom
private:
static double GammaFctHalfInt( double x);
constexpr static unsigned int factorial( unsigned int n);
const int nu; // degrees of freedom
};
inline ChiSquared::ChiSquared( const int nu ) : nu(nu) {
assert( nu>=0 && "degrees of freedom negative");
}
/*!
* \brief probability density function
* \param x real number [0,inf)
* \return value of the probability density function (pdf)
*/
inline double ChiSquared::pdf( const double x) const
{
if ( x<0 ) {
return 0.;
}
// K.-R. Koch (261.1)
return 1. / ( pow(2,0.5*nu) * GammaFctHalfInt(0.5*nu) )
* pow(x, 0.5*nu-1) * exp(-0.5*x);
//qDebug().noquote() << QString("chi2pdf(%1,%2) = %3").arg(x).arg(nu).arg(f);
}
/*!
* \brief cumulative density function
* \param x real number
* \return value of cumulative density function
*/
inline Prob ChiSquared::cdf( const double x) const
{
if ( x<0 ) {
return Prob(0);
}
// K.-R. Koch, (261.6).
constexpr int num_iter_max = 20;
constexpr double threshold = 1e-6;
double sum = 0.;
double prod = 1.;
for (int i=1; i<num_iter_max; i++) {
prod *= static_cast<double>( nu+2*i );
const double summand = pow(x, i) / prod;
sum += summand;
if ( summand < threshold ) {
break;
}
}
return Prob(
pow( 0.5*x, 0.5*nu) * exp(-0.5*x) / GammaFctHalfInt( 0.5*(nu+2) )*(1. +sum )
);
//qDebug().noquote() << QString("chi2cdf(%1, %2) = %3").arg(x).arg(nu).arg(F);
}
/*!
* \brief Quantile of the chi-squared distribution.
* \param P probability [0,1]
* \return Quantile x so that cdf(x,nu) = P.
*/
inline double ChiSquared::icdf( const Prob P) const
{
const double alpha = 1.-P();
// approx., K.-R. Koch (261.11)
// const StandardNormal normal;
const double xa = Stats::StandardNormal::icdf(P);
const auto df = static_cast<double>(nu);
const double t = 2./(9.*df);
double q = df* pow( xa*sqrt( t )+1.-t, 3 );
// iterative refinement of q, K.-R. Koch (261.13)
constexpr int numIterMax = 10;
constexpr double threshold = 1e-7;
for ( int i=0; i<numIterMax; i++) {
const double alpha_n = 1. - cdf(q)();
const double summand = (alpha_n - alpha) / pdf(q);
q += summand;
if ( fabs(summand) < threshold ) {
break;
}
}
//qDebug().noquote() << QString("chi2inv(%1,%2) = %3").arg(P).arg(nu).arg(q);
return q;
}
inline double ChiSquared::rnd() const
{
const Stats::Gamma gamma( nu/2., 2.);
return gamma.rnd();
}
/*!
* \brief Gamma function for arguments being multiples of 0.5
* \param x parameter
* \return GammaFct(x)
*/
inline double ChiSquared::GammaFctHalfInt( const double x)
{
assert( trunc(2*x)==2*x && "argument not a multiple of 1/2" );
// K.-R. Koch, (243.5), (243.6).
constexpr double sqrt_pi = 1.772453850905516; // = sqrt(3.14159);
// case 1, 2, 3, ...
if ( trunc(x)==x ) {
return factorial( static_cast<unsigned int>(x)-1); // x \in N, x>0
}
// else: case 0.5, 1.5, 2.5, ...
const int p = static_cast<int>( x-0.5 );
int pr = 1; // product
for ( int i=1; i<=2*p-1; i+=2 ) {
pr = pr*i;
}
return pr*sqrt_pi/static_cast<double>( pow(2,p) );
}
/*!
* \brief factorial, compile-time enabled evaluation
* \param n non-negative integer [0,17]
* \return factorial n!
*/
constexpr unsigned int ChiSquared::factorial( const unsigned int n)
{
// prevent silliness and overflow:
assert( n<18 ); // 17! = 3.5569e+14
unsigned int x = 1; // 0! = 1
for ( unsigned int i=2; i<=n; ++i) {
x *= i;
}
return x;
}
} // namespace Stats
#endif // CHISQUARED_H