mir.stat.distribution.gev source code

1 /++
2 This module contains algorithms for the $(LINK2 https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution, Generalized Extreme Value (GEV) 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 
12 module mir.stat.distribution.gev;
13 
14 import mir.internal.utility: isFloatingPoint;
15 
16 import mir.math.common: fabs, exp, pow, log;
17 
18 /++
19 Computes the generalized extreme value (GEV) probability density function (PDF).
20 
21 Params:
22     x = value to evaluate
23     mu = location
24     sigma = scale
25     xi = shape
26 
27 See_also:
28     $(LINK2 https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution, Generalized Extreme Value (GEV) Distribution)
29 +/
30 @safe pure nothrow @nogc
31 T gevPDF(T)(const T x, const T mu, const T sigma, const T xi)
32     if (isFloatingPoint!T)
33     in (xi >= 0 || x <= mu - sigma / xi, "if xi is less than zero, x must be less than or equal to mu - sigma / xi")
34     in (xi <= 0 || x >= mu - sigma / xi, "if xi is greater than zero, xi must be greater than or equal to mu - sigma / xi")
35 {
36     auto s = (x - mu) / sigma;
37     if (xi.fabs <= T.min_normal)
38     {
39         auto t = exp(-s);
40         return t * exp(-t);
41     }
42     auto v = 1 + xi * s;
43     if (v <= 0)
44         return 0;
45     auto a = pow(v, -1 / xi);
46     return a * exp(-a) / (v * sigma);
47 }
48 
49 ///
50 version(mir_stat_test)
51 @safe pure nothrow @nogc
52 unittest
53 {
54     import mir.test: shouldApprox;
55 
56     gevPDF(-3, 2, 3, -0.5).shouldApprox == 0.02120353011709564;
57     gevPDF(-1, 2, 3, +0.5).shouldApprox == 0.04884170370329114;
58     gevPDF(-1, 2, 3, 0.0).shouldApprox == 0.1793740787340172;
59 }
60 
61 // Checking v <= 0 branch
62 version(mir_stat_test)
63 @safe pure nothrow @nogc
64 unittest
65 {
66     import mir.test: should;
67     gevPDF(-1.0, 0, 1, 1).should == 0;
68 }
69 
70 /++
71 Computes the generalized extreme value (GEV) cumulatve distribution function (CDF).
72 
73 Params:
74     x = value to evaluate
75     mu = location
76     sigma = scale
77     xi = shape
78 
79 See_also:
80     $(LINK2 https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution, Generalized Extreme Value (GEV) Distribution)
81 +/
82 @safe pure nothrow @nogc
83 T gevCDF(T)(const T x, const T mu, const T sigma, const T xi)
84     if (isFloatingPoint!T)
85     in (xi >= 0 || x <= mu - sigma / xi, "if xi is less than zero, x must be less than or equal to mu - sigma / xi")
86     in (xi <= 0 || x >= mu - sigma / xi, "if xi is greater than zero, xi must be greater than or equal to mu - sigma / xi")
87 {
88     auto s = (x - mu) / sigma;
89     if (xi.fabs <= T.min_normal)
90         return exp(-exp(-s));
91     auto v = 1 + xi * s;
92     if (v <= 0)
93         return xi > 0 ? 0 : 1;
94     auto a = pow(v, -1 / xi);
95     return exp(-a);
96 }
97 
98 ///
99 version(mir_stat_test)
100 @safe pure nothrow @nogc
101 unittest
102 {
103     import mir.test: shouldApprox;
104 
105     gevCDF(-3, 2, 3, -0.5).shouldApprox == 0.034696685646156494;
106     gevCDF(-1, 2, 3, +0.5).shouldApprox == 0.01831563888873418;
107     gevCDF(-1, 2, 3, 0.0).shouldApprox == 0.06598803584531254;
108 }
109 
110 // Checking v <= 0 branch
111 version(mir_stat_test)
112 @safe pure nothrow @nogc
113 unittest
114 {
115     import mir.test: should;
116     gevCDF(-1.0, 0, 1, 1).should == 0;
117     gevCDF(1.0, 0, 1, -1).should == 1;
118 }
119 
120 /++
121 Computes the generalized extreme value (GEV) complementary cumulatve distribution function (CCDF).
122 
123 Params:
124     x = value to evaluate
125     mu = location
126     sigma = scale
127     xi = shape
128 
129 See_also:
130     $(LINK2 https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution, Generalized Extreme Value (GEV) Distribution)
131 +/
132 @safe pure nothrow @nogc
133 T gevCCDF(T)(const T x, const T mu, const T sigma, const T xi)
134     if (isFloatingPoint!T)
135     in (xi >= 0 || x <= mu - sigma / xi, "if xi is less than zero, x must be less than or equal to mu - sigma / xi")
136     in (xi <= 0 || x >= mu - sigma / xi, "if xi is greater than zero, xi must be greater than or equal to mu - sigma / xi")
137 {
138     return 1 - gevCDF(x, mu, sigma, xi);
139 }
140 
141 ///
142 version(mir_stat_test)
143 @safe pure nothrow @nogc
144 unittest
145 {
146     import mir.test: shouldApprox;
147 
148     gevCCDF(-3, 2, 3, -0.5).shouldApprox == 0.965303314353844;
149     gevCCDF(-1, 2, 3, +0.5).shouldApprox == 0.981684361111266;
150     gevCCDF(-1, 2, 3, 0.0).shouldApprox == 0.934011964154687;
151 }
152 
153 // Checking v <= 0 branch
154 version(mir_stat_test)
155 @safe pure nothrow @nogc
156 unittest
157 {
158     import mir.test: should;
159     gevCCDF(-1.0, 0, 1, 1).should == 1;
160     gevCCDF(1.0, 0, 1, -1).should == 0;
161 }
162 
163 /++
164 Computes the generalized extreme value (GEV) inverse cumulative distribution function (InvCDF).
165 
166 Params:
167     p = value to evaluate
168     mu = location
169     sigma = scale
170     xi = shape
171 
172 See_also:
173     $(LINK2 https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution, Generalized Extreme Value (GEV) Distribution)
174 +/
175 @safe pure nothrow @nogc
176 T gevInvCDF(T)(const T p, const T mu, const T sigma, const T xi)
177     if (isFloatingPoint!T)
178     in (p >= 0, "p must be greater than or equal to 0")
179     in (p <= 1, "p must be less than or equal to 1")
180 {
181     auto logp = log(p);
182     if (xi.fabs <= T.min_normal)
183         return mu - sigma * log(-logp);
184     return mu + (pow(-logp, -xi) - 1) * sigma / xi;
185 }
186 
187 ///
188 version(mir_stat_test)
189 @safe pure nothrow @nogc
190 unittest
191 {
192     import mir.test: shouldApprox;
193 
194     gevInvCDF(0.034696685646156494, 2, 3, -0.5).shouldApprox == -3;
195     gevInvCDF(0.01831563888873418, 2, 3, +0.5).shouldApprox == -1;
196     gevInvCDF(0.06598803584531254, 2, 3, 0.0).shouldApprox == -1;
197 }
198 
199 /++
200 Computes the generalized extreme value (GEV) log probability density function (LPDF).
201 
202 Params:
203     x = value to evaluate
204     mu = location
205     sigma = scale
206     xi = shape
207 
208 See_also:
209     $(LINK2 https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution, Generalized Extreme Value (GEV) Distribution)
210 +/
211 @safe pure nothrow @nogc
212 T gevLPDF(T)(const T x, const T mu, const T sigma, const T xi)
213     if (isFloatingPoint!T)
214     in (xi >= 0 || x <= mu - sigma / xi, "if xi is less than zero, x must be less than or equal to mu - sigma / xi")
215     in (xi <= 0 || x >= mu - sigma / xi, "if xi is greater than zero, xi must be greater than or equal to mu - sigma / xi")
216 {
217     import mir.math.common: log;
218 
219     auto s = (x - mu) / sigma;
220     if (xi.fabs <= T.min_normal)
221     {
222         auto t = exp(-s);
223         return log(t) - t;
224     }
225     auto v = 1 + xi * s;
226     if (v <= 0)
227         return -double.infinity;
228     auto a = pow(v, -1 / xi);
229     return log(a) - a - log(v * sigma);
230 }
231 
232 ///
233 version(mir_stat_test)
234 @safe pure nothrow @nogc
235 unittest
236 {
237     import mir.test: shouldApprox;
238 
239     gevLPDF(-3, 2, 3, -0.5).shouldApprox == -3.85358759620891;
240     gevLPDF(-1, 2, 3, +0.5).shouldApprox == -3.01917074698827;
241     gevLPDF(-1, 2, 3, 0.0).shouldApprox == -1.71828182845905;
242 }
243 
244 // Checking v <= 0 branch
245 version(mir_stat_test)
246 @safe pure nothrow @nogc
247 unittest
248 {
249     import mir.test: shouldApprox;
250     gevLPDF(-1.0, 0, 1, 1).shouldApprox == -double.infinity;
251 }
The OpenD Programming Language
