mir.stat.distribution.chi2 source code

1 /++
2 This module contains algorithms for the $(LINK2 https://en.wikipedia.org/wiki/Chi-squared_distribution, Chi-squared Distribution).
3 
4 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
5 
6 Authors: Ilia Ki, John Michael Hall
7 
8 Copyright: 2022-3 Mir Stat Authors.
9 
10 +/
11 module mir.stat.distribution.chi2;
12 
13 import mir.internal.utility: isFloatingPoint;
14 
15 /++
16 Computes the Chi-squared probability density function (PDF).
17 
18 Params:
19     x = value to evaluate PDF
20     k = degrees of freedom
21 
22 See_also:
23     $(LINK2 https://en.wikipedia.org/wiki/Chi-squared_distribution, Chi-squared Distribution)
24 +/
25 @safe pure nothrow @nogc
26 T chi2PDF(T)(const T x, const uint k)
27     if (isFloatingPoint!T)
28     in (x >= 0, "x must be greater than or equal to 0")
29     in (k >= 1, "k must be greater than or equal to 1")
30 {
31     import mir.stat.distribution.gamma: gammaPDF;
32 
33     return gammaPDF(x, T(k) * 0.5f, 2);
34 }
35 
36 ///
37 version(mir_stat_test)
38 @safe pure nothrow @nogc
39 unittest
40 {
41     import mir.test: shouldApprox;
42     0.2.chi2PDF(2).shouldApprox == 0.4524187;
43 }
44 
45 /++
46 Computes the Chi-squared cumulative distribution function (CDF).
47 
48 Params:
49     x = value to evaluate CDF
50     k = degrees of freedom
51 
52 See_also:
53     $(LINK2 https://en.wikipedia.org/wiki/Chi-squared_distribution, Chi-squared Distribution)
54 +/
55 @safe pure nothrow @nogc
56 T chi2CDF(T)(const T x, const uint k)
57     if (isFloatingPoint!T)
58     in (x >= 0, "x must be greater than or equal to 0")
59     in (k >= 1, "k must be greater than or equal to 1")
60 {
61     import mir.stat.distribution.gamma: gammaCDF;
62 
63     return gammaCDF(x, T(k) * 0.5f, 2);
64 }
65 
66 ///
67 version(mir_stat_test)
68 @safe pure nothrow @nogc
69 unittest
70 {
71     import mir.test: shouldApprox;
72     0.2.chi2CDF(2).shouldApprox == 0.09516258;
73 }
74 
75 /++
76 Computes the Chi-squared complementary cumulative distribution function (CCDF).
77 
78 Params:
79     x = value to evaluate CCDF
80     k = degrees of freedom
81 
82 See_also:
83     $(LINK2 https://en.wikipedia.org/wiki/Chi-squared_distribution, Chi-squared Distribution)
84 +/
85 @safe pure nothrow @nogc
86 T chi2CCDF(T)(const T x, const uint k)
87     if (isFloatingPoint!T)
88     in (x >= 0, "x must be greater than or equal to 0")
89     in (k >= 1, "k must be greater than or equal to 1")
90 {
91     import mir.stat.distribution.gamma: gammaCCDF;
92 
93     return gammaCCDF(x, T(k) * 0.5f , 2);
94 }
95 
96 ///
97 version(mir_stat_test)
98 @safe pure nothrow @nogc
99 unittest
100 {
101     import mir.test: shouldApprox;
102     0.2.chi2CCDF(2).shouldApprox == 0.9048374;
103 }
104 
105 /++
106 Computes the Chi-squared inverse cumulative distribution function (InvCDF).
107 
108 Params:
109     p = value to evaluate InvCDF
110     k = degrees of freedom
111 
112 See_also:
113     $(LINK2 https://en.wikipedia.org/wiki/Chi-squared_distribution, Chi-squared Distribution)
114 +/
115 @safe pure nothrow @nogc
116 T chi2InvCDF(T)(const T p, const uint k)
117     if (isFloatingPoint!T)
118     in (p >= 0, "p must be greater than or equal to 0")
119     in (p <= 1, "p must be less than or equal to 1")
120     in (k >= 1, "k must be greater than or equal to 1")
121 {
122     import mir.stat.distribution.gamma: gammaInvCDF;
123 
124     return gammaInvCDF(p, T(k) * 0.5f, 2);
125 }
126 
127 ///
128 version(mir_stat_test)
129 @safe pure nothrow @nogc
130 unittest
131 {
132     import mir.test: shouldApprox;
133     0.09516258.chi2InvCDF(2).shouldApprox == 0.2;
134 }
135 
136 /++
137 Computes the Chi-squared probability density function (LPDF).
138 
139 Params:
140     x = value to evaluate LPDF
141     k = degrees of freedom
142 
143 See_also:
144     $(LINK2 https://en.wikipedia.org/wiki/Chi-squared_distribution, Chi-squared Distribution)
145 +/
146 @safe pure nothrow @nogc
147 T chi2LPDF(T)(const T x, const uint k)
148     if (isFloatingPoint!T)
149     in (x >= 0, "x must be greater than or equal to 0")
150     in (k >= 1, "k must be greater than or equal to 1")
151 {
152     import mir.math.common: log;
153     import mir.math.constant: LN2;
154     import std.mathspecial: logGamma;
155 
156     if (x == 0) {
157         if (k > 2) {
158             return -T.infinity;
159         } else if (k < 2) {
160             return T.infinity;
161         } else {
162             return -T(LN2);
163         }
164     }
165 
166     return (T(k) * 0.5f - 1) * (log(x) - T(LN2)) - x *0.5f - cast(T) logGamma(T(k) * 0.5f) - T(LN2);
167 }
168 
169 ///
170 version(mir_stat_test)
171 @safe pure nothrow @nogc
172 unittest
173 {
174     import mir.test: shouldApprox;
175     0.2.chi2LPDF(2).shouldApprox == -0.7931472;
176 }
177 
178 //
179 version(mir_stat_test)
180 @safe pure nothrow @nogc
181 unittest
182 {
183     import mir.test: shouldApprox;
184     0.0.chi2LPDF(2).shouldApprox == -0.6931472;
185     0.0.chi2LPDF(1).shouldApprox == double.infinity;
186     0.0.chi2LPDF(3).shouldApprox == -double.infinity;
187 }
The OpenD Programming Language
