mir.stat.distribution.cauchy source code

1 /++
2 This module contains algorithms for the $(LINK2 https://en.wikipedia.org/wiki/Cauchy_distribution, Cauchy Distribution).
3 
4 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
5 
6 Authors: John Michael Hall
7 
8 Copyright: 2023 Mir Stat Authors.
9 
10 +/
11 
12 module mir.stat.distribution.cauchy;
13 
14 import mir.internal.utility: isFloatingPoint;
15 
16 /++
17 Computes the Cauchy probability density function (PDF).
18 
19 Params:
20     x = value to evaluate PDF
21 
22 See_also:
23     $(LINK2 https://en.wikipedia.org/wiki/Cauchy_distribution, Cauchy Distribution)
24 +/
25 @safe pure nothrow @nogc
26 T cauchyPDF(T)(const T x)
27     if (isFloatingPoint!T)
28 {
29     import mir.math.constant: M_1_PI;
30 
31     return T(M_1_PI) / (1 + x * x);
32 }
33 
34 /++
35 Ditto, with location and scale parameters (by standardizing `x`).
36 
37 Params:
38     x = value to evaluate PDF
39     location = location parameter
40     scale = scale parameter
41 +/
42 @safe pure nothrow @nogc
43 T cauchyPDF(T)(const T x, const T location, const T scale)
44     if (isFloatingPoint!T)
45     in (scale > 0, "scale must be greater than zero")
46 {
47     return cauchyPDF((x - location) / scale) / scale;
48 }
49 
50 ///
51 version(mir_stat_test)
52 @safe pure nothrow @nogc
53 unittest {
54     import mir.test: shouldApprox;
55 
56     cauchyPDF(-3.0).shouldApprox == 0.03183099;
57     cauchyPDF(-2.0).shouldApprox == 0.06366198;
58     cauchyPDF(-1.0).shouldApprox == 0.1591549;
59     cauchyPDF(0.0).shouldApprox == 0.3183099;
60     cauchyPDF(1.0).shouldApprox == 0.1591549;
61     cauchyPDF(2.0).shouldApprox == 0.06366198;
62     cauchyPDF(3.0).shouldApprox == 0.03183099;
63 
64     // Can include location/scale
65     cauchyPDF(-3.0, 1, 2).shouldApprox == 0.03183099;
66     cauchyPDF(-2.0, 1, 2).shouldApprox == 0.04897075;
67     cauchyPDF(-1.0, 1, 2).shouldApprox == 0.07957747;
68     cauchyPDF(0.0, 1, 2).shouldApprox == 0.127324;
69     cauchyPDF(1.0, 1, 2).shouldApprox == 0.1591549;
70     cauchyPDF(2.0, 1, 2).shouldApprox == 0.127324;
71     cauchyPDF(3.0, 1, 2).shouldApprox == 0.07957747;
72 }
73 
74 /++
75 Computes the Cauchy cumulative distribution function (CDF).
76 
77 Params:
78     x = value to evaluate CDF
79 
80 See_also:
81     $(LINK2 https://en.wikipedia.org/wiki/Cauchy_distribution, Cauchy Distribution)
82 +/
83 @safe pure nothrow @nogc
84 T cauchyCDF(T)(const T x)
85     if (isFloatingPoint!T)
86 {
87     import mir.math.constant: M_1_PI;
88     import std.math.trigonometry: atan;
89 
90     return 0.5 + T(M_1_PI) * atan(x);
91 }
92 
93 /++
94 Ditto, with location and scale parameters (by standardizing `x`).
95 
96 Params:
97     x = value to evaluate CDF
98     location = location parameter
99     scale = scale parameter
100 +/
101 @safe pure nothrow @nogc
102 T cauchyCDF(T)(const T x, const T location, const T scale)
103     if (isFloatingPoint!T)
104     in (scale > 0, "scale must be greater than zero")
105 {
106     return cauchyCDF((x - location) / scale);
107 }
108 
109 ///
110 version(mir_stat_test)
111 @safe pure nothrow @nogc
112 unittest {
113     import mir.test: shouldApprox;
114 
115     cauchyCDF(-3.0).shouldApprox == 0.1024164;
116     cauchyCDF(-2.0).shouldApprox == 0.1475836;
117     cauchyCDF(-1.0).shouldApprox == 0.25;
118     cauchyCDF(0.0).shouldApprox == 0.5;
119     cauchyCDF(1.0).shouldApprox == 0.75;
120     cauchyCDF(2.0).shouldApprox == 0.8524164;
121     cauchyCDF(3.0).shouldApprox == 0.8975836;
122 
123     // Can include location/scale
124     cauchyCDF(-3.0, 1, 2).shouldApprox == 0.1475836;
125     cauchyCDF(-2.0, 1, 2).shouldApprox == 0.187167;
126     cauchyCDF(-1.0, 1, 2).shouldApprox == 0.25;
127     cauchyCDF(0.0, 1, 2).shouldApprox == 0.3524164;
128     cauchyCDF(1.0, 1, 2).shouldApprox == 0.5;
129     cauchyCDF(2.0, 1, 2).shouldApprox == 0.6475836;
130     cauchyCDF(3.0, 1, 2).shouldApprox == 0.75;
131 }
132 
133 /++
134 Computes the Cauchy complementary cumulative distribution function (CCDF).
135 
136 Params:
137     x = value to evaluate CCDF
138 
139 See_also:
140     $(LINK2 https://en.wikipedia.org/wiki/Cauchy_distribution, Cauchy Distribution)
141 +/
142 @safe pure nothrow @nogc
143 T cauchyCCDF(T)(const T x)
144     if (isFloatingPoint!T)
145 {
146     return cauchyCDF(-x);
147 }
148 
149 /++
150 Ditto, with location and scale parameters (by standardizing `x`).
151 
152 Params:
153     x = value to evaluate CCDF
154     location = location parameter
155     scale = scale parameter
156 +/
157 @safe pure nothrow @nogc
158 T cauchyCCDF(T)(const T x, const T location, const T scale)
159     if (isFloatingPoint!T)
160     in (scale > 0, "scale must be greater than zero")
161 {
162     return cauchyCDF((location - x) / scale);
163 }
164 
165 ///
166 version(mir_stat_test)
167 @safe pure nothrow @nogc
168 unittest {
169     import mir.test: shouldApprox;
170 
171     cauchyCCDF(-3.0).shouldApprox == 0.8975836;
172     cauchyCCDF(-2.0).shouldApprox == 0.8524164;
173     cauchyCCDF(-1.0).shouldApprox == 0.75;
174     cauchyCCDF(0.0).shouldApprox == 0.5;
175     cauchyCCDF(1.0).shouldApprox == 0.25;
176     cauchyCCDF(2.0).shouldApprox == 0.1475836;
177     cauchyCCDF(3.0).shouldApprox == 0.1024164;
178 
179     // Can include location/scale
180     cauchyCCDF(-3.0, 1, 2).shouldApprox == 0.8524164;
181     cauchyCCDF(-2.0, 1, 2).shouldApprox == 0.812833;
182     cauchyCCDF(-1.0, 1, 2).shouldApprox == 0.75;
183     cauchyCCDF(0.0, 1, 2).shouldApprox == 0.6475836;
184     cauchyCCDF(1.0, 1, 2).shouldApprox == 0.5;
185     cauchyCCDF(2.0, 1, 2).shouldApprox == 0.3524164;
186     cauchyCCDF(3.0, 1, 2).shouldApprox == 0.25;
187 }
188 
189 /++
190 Computes the Cauchy inverse cumulative distribution function (InvCDF).
191 
192 Params:
193     p = value to evaluate InvCDF
194 
195 See_also:
196     $(LINK2 https://en.wikipedia.org/wiki/Cauchy_distribution, Cauchy Distribution)
197 +/
198 @safe pure nothrow @nogc
199 T cauchyInvCDF(T)(const T p)
200     if (isFloatingPoint!T)
201     in (p >= 0, "p must be greater than or equal to 0")
202     in (p <= 1, "p must be less than or equal to 1")
203 {
204     import mir.math.constant: PI;
205     import std.math.trigonometry: tan;
206 
207     if (p > 0 && p < 1) {
208         return tan(T(PI) * (p - 0.5));
209     } else if (p == 0) {
210         return -T.infinity;
211     } else if (p == 1) {
212         return T.infinity;
213     }
214     assert(0, "Should not be here");
215 }
216 
217 /++
218 Ditto, with location and scale parameters.
219 
220 Params:
221     p = value to evaluate InvCDF
222     location = location parameter
223     scale = scale parameter
224 +/
225 @safe pure nothrow @nogc
226 T cauchyInvCDF(T)(const T p, const T location, const T scale)
227     if (isFloatingPoint!T)
228     in (p >= 0, "p must be greater than or equal to 0")
229     in (p <= 1, "p must be less than or equal to 1")
230     in (scale > 0, "scale must be greater than zero")
231 {
232      return location + scale * cauchyInvCDF(p);
233 }
234 
235 ///
236 version(mir_stat_test)
237 @safe pure nothrow @nogc
238 unittest {
239     import mir.test: shouldApprox;
240 
241     cauchyInvCDF(0.0).shouldApprox == -double.infinity;
242     cauchyInvCDF(0.1).shouldApprox == -3.077684;
243     cauchyInvCDF(0.2).shouldApprox == -1.376382;
244     cauchyInvCDF(0.3).shouldApprox == -0.7265425;
245     cauchyInvCDF(0.4).shouldApprox == -0.3249197;
246     cauchyInvCDF(0.5).shouldApprox == 0.0;
247     cauchyInvCDF(0.6).shouldApprox == 0.3249197;
248     cauchyInvCDF(0.7).shouldApprox == 0.7265425;
249     cauchyInvCDF(0.8).shouldApprox == 1.376382;
250     cauchyInvCDF(0.9).shouldApprox == 3.077684;
251     cauchyInvCDF(1.0).shouldApprox == double.infinity;
252 
253     // Can include location/scale
254     cauchyInvCDF(0.2, 1, 2).shouldApprox == -1.752764;
255     cauchyInvCDF(0.4, 1, 2).shouldApprox == 0.3501606;
256     cauchyInvCDF(0.6, 1, 2).shouldApprox == 1.649839;
257     cauchyInvCDF(0.8, 1, 2).shouldApprox == 3.752764;
258 }
259 
260 /++
261 Computes the Cauchy log probability density function (LPDF).
262 
263 Params:
264     x = value to evaluate LPDF
265 
266 See_also:
267     $(LINK2 https://en.wikipedia.org/wiki/Cauchy_distribution, Cauchy Distribution)
268 +/
269 @safe pure nothrow @nogc
270 T cauchyLPDF(T)(const T x)
271     if (isFloatingPoint!T)
272 {
273     import mir.math.common: log;
274     import mir.stat.constant: LOGPI;
275 
276     return -T(LOGPI) - log(1 + x * x);
277 }
278 
279 /++
280 Ditto, with location and scale parameters (by standardizing `x`).
281 
282 Params:
283     x = value to evaluate LPDF
284     location = location parameter
285     scale = scale parameter
286 +/
287 @safe pure nothrow @nogc
288 T cauchyLPDF(T)(const T x, const T location, const T scale)
289     if (isFloatingPoint!T)
290     in (scale > 0, "scale must be greater than zero")
291 {
292     import mir.math.common: log;
293 
294     return cauchyLPDF((x - location) / scale) - log(scale);
295 }
296 
297 ///
298 version(mir_stat_test)
299 @safe pure nothrow @nogc
300 unittest {
301     import mir.math.common: log;
302     import mir.test: shouldApprox;
303 
304     cauchyLPDF(-3.0).shouldApprox == log(0.03183099);
305     cauchyLPDF(-2.0).shouldApprox == log(0.06366198);
306     cauchyLPDF(-1.0).shouldApprox == log(0.1591549);
307     cauchyLPDF(0.0).shouldApprox == log(0.3183099);
308     cauchyLPDF(1.0).shouldApprox == log(0.1591549);
309     cauchyLPDF(2.0).shouldApprox == log(0.06366198);
310     cauchyLPDF(3.0).shouldApprox == log(0.03183099);
311 
312     // Can include location/scale
313     cauchyLPDF(-3.0, 1, 2).shouldApprox == log(0.03183099);
314     cauchyLPDF(-2.0, 1, 2).shouldApprox == log(0.04897075);
315     cauchyLPDF(-1.0, 1, 2).shouldApprox == log(0.07957747);
316     cauchyLPDF(0.0, 1, 2).shouldApprox == log(0.127324);
317     cauchyLPDF(1.0, 1, 2).shouldApprox == log(0.1591549);
318     cauchyLPDF(2.0, 1, 2).shouldApprox == log(0.127324);
319     cauchyLPDF(3.0, 1, 2).shouldApprox == log(0.07957747);
320 }
The OpenD Programming Language
