mir.stat.distribution.generalized

1 /++
2 This module contains algorithms for the $(LINK2 https://en.wikipedia.org/wiki/Generalized_Pareto_distribution, Generalized Pareto Distribution).
3 
4 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
5 
6 Authors: John Michael Hall
7 
8 Copyright: 2022-3 Mir Stat Authors.
9 
10 +/
11 
12 module mir.stat.distribution.generalized_pareto;
13 
14 import mir.internal.utility: isFloatingPoint;
15 
16 /++
17 Computes the generalized pareto probability density function (PDF).
18 
19 Params:
20     x = value to evaluate PDF
21     mu = location parameter
22     sigma = scale parameter
23     xi = shape parameter
24 
25 See_also:
26     $(LINK2 https://en.wikipedia.org/wiki/Generalized_Pareto_distribution, Generalized Pareto Distribution)
27 +/
28 @safe pure nothrow @nogc
29 T generalizedParetoPDF(T)(const T x, const T mu, const T sigma, const T xi)
30     if (isFloatingPoint!T)
31     in (sigma > 0, "sigma must be greater than zero")
32     in (x >= mu, "x must be greater than or equal to mu")
33     in (xi >= 0 || (xi < 0 && x <= (mu - sigma / xi)), "if xi is less than zero, x must be less than mu - sigma / xi")
34 {
35     import mir.math.common: exp, pow;
36 
37     const T z = (x - mu) / sigma;
38     if (xi != 0) {
39         return (cast(T) 1 / sigma) * pow(1 + xi * z, -(cast(T) 1 / xi + 1));
40     } else {
41         return exp(-z);
42     }
43 }
44 
45 ///
46 version(mir_stat_test)
47 @safe pure nothrow @nogc
48 unittest {
49     import mir.test: shouldApprox;
50 
51     1.0.generalizedParetoPDF(1, 1, 0.5).shouldApprox == 1;
52     2.0.generalizedParetoPDF(1, 1, 0.5).shouldApprox == 0.2962963;
53     3.0.generalizedParetoPDF(2, 3, 0.25).shouldApprox == 0.2233923;
54     5.0.generalizedParetoPDF(2, 3, 0).shouldApprox == 0.3678794;
55 }
56 
57 /++
58 Computes the generalized pareto cumulative distribution function (CDF).
59 
60 Params:
61     x = value to evaluate CDF
62     mu = location parameter
63     sigma = scale parameter
64     xi = shape parameter
65 
66 See_also:
67     $(LINK2 https://en.wikipedia.org/wiki/Generalized_Pareto_distribution, Generalized Pareto Distribution)
68 +/
69 @safe pure nothrow @nogc
70 T generalizedParetoCDF(T)(const T x, const T mu, const T sigma, const T xi)
71     if (isFloatingPoint!T)
72     in (sigma > 0, "sigma must be greater than zero")
73     in (x >= mu, "x must be greater than or equal to mu")
74     in (xi >= 0 || (xi < 0 && x <= (mu - sigma / xi)), "if xi is less than zero, x must be less than mu - sigma / xi")
75 {
76     import mir.math.common: exp, pow;
77 
78     const T z = (x - mu) / sigma;
79     if (xi != 0) {
80         return 1 - pow(1 + xi * z, -(cast(T) 1) / xi);
81     } else {
82         return 1 - exp(-z);
83     }
84 }
85 
86 ///
87 version(mir_stat_test)
88 @safe pure nothrow @nogc
89 unittest {
90     import mir.test: shouldApprox;
91 
92     1.0.generalizedParetoCDF(1, 1, 0.5).shouldApprox == 0;
93     2.0.generalizedParetoCDF(1, 1, 0.5).shouldApprox == 0.5555556;
94     3.0.generalizedParetoCDF(2, 3, 0.25).shouldApprox == 0.273975;
95     5.0.generalizedParetoCDF(2, 3, 0).shouldApprox == 0.6321206;
96 }
97 
98 /++
99 Computes the generalized pareto complementary cumulative distribution function (CCDF).
100 
101 Params:
102     x = value to evaluate CCDF
103     mu = location parameter
104     sigma = scale parameter
105     xi = shape parameter
106 
107 See_also:
108     $(LINK2 https://en.wikipedia.org/wiki/Generalized_Pareto_distribution, Generalized Pareto Distribution)
109 +/
110 @safe pure nothrow @nogc
111 T generalizedParetoCCDF(T)(const T x, const T mu, const T sigma, const T xi)
112     if (isFloatingPoint!T)
113     in (sigma > 0, "sigma must be greater than zero")
114     in (x >= mu, "x must be greater than or equal to mu")
115     in (xi >= 0 || (xi < 0 && x <= (mu - sigma / xi)), "if xi is less than zero, x must be less than mu - sigma / xi")
116 {
117     import mir.math.common: exp, pow;
118 
119     const T z = (x - mu) / sigma;
120     if (xi != 0) {
121         return pow(1 + xi * z, -(cast(T) 1) / xi);
122     } else {
123         return exp(-z);
124     }
125 }
126 
127 ///
128 version(mir_stat_test)
129 @safe pure nothrow @nogc
130 unittest {
131     import mir.test: shouldApprox;
132 
133     1.0.generalizedParetoCCDF(1, 1, 0.5).shouldApprox == 1;
134     2.0.generalizedParetoCCDF(1, 1, 0.5).shouldApprox == 0.4444444;
135     3.0.generalizedParetoCCDF(2, 3, 0.25).shouldApprox == 0.726025;
136     5.0.generalizedParetoCCDF(2, 3, 0).shouldApprox == 0.3678794;
137 }
138 
139 /++
140 Computes the generalized pareto inverse cumulative distribution function (InvCDF).
141 
142 Params:
143     p = value to evaluate InvCDF
144     mu = location parameter
145     sigma = scale parameter
146     xi = shape parameter
147 
148 See_also:
149     $(LINK2 https://en.wikipedia.org/wiki/Generalized_Pareto_distribution, Generalized Pareto Distribution)
150 +/
151 @safe pure nothrow @nogc
152 T generalizedParetoInvCDF(T)(const T p, const T mu, const T sigma, const T xi)
153     if (isFloatingPoint!T)
154     in (p >= 0, "p must be greater than or equal to 0")
155     in (p <= 1, "p must be less than or equal to 1")
156     in (sigma > 0, "sigma must be greater than zero")
157 {
158     import mir.math.common: pow;
159     import mir.math.internal.log1p: log1p;
160 
161     T output;
162     if (xi != 0) {
163         output = (cast(T) 1 / xi) * (pow(1 - p, -xi) - 1);
164     } else {
165         output = -log1p(-p);
166     }
167     return mu + sigma * output;
168 }
169 
170 ///
171 version(mir_stat_test)
172 @safe pure nothrow @nogc
173 unittest {
174     import mir.test: shouldApprox;
175 
176     0.0.generalizedParetoInvCDF(1, 1, 0.5).shouldApprox == 1;
177     0.5555556.generalizedParetoInvCDF(1, 1, 0.5).shouldApprox == 2;
178     0.273975.generalizedParetoInvCDF(2, 3, 0.25).shouldApprox == 3;
179     0.6321206.generalizedParetoInvCDF(2, 3, 0).shouldApprox == 5;    
180 }
181 
182 /++
183 Computes the generalized pareto log probability density function (LPDF).
184 
185 Params:
186     x = value to evaluate LPDF
187     mu = location parameter
188     sigma = scale parameter
189     xi = shape parameter
190 
191 See_also:
192     $(LINK2 https://en.wikipedia.org/wiki/Generalized_Pareto_distribution, Generalized Pareto Distribution)
193 +/
194 @safe pure nothrow @nogc
195 T generalizedParetoLPDF(T)(const T x, const T mu, const T sigma, const T xi)
196     if (isFloatingPoint!T)
197     in (sigma > 0, "sigma must be greater than zero")
198     in (x >= mu, "x must be greater than or equal to mu")
199     in (xi >= 0 || (xi < 0 && x <= (mu - sigma / xi)), "if xi is less than zero, x must be less than mu - sigma / xi")
200 {
201     import mir.math.common: log;
202     import mir.math.internal.xlogy: xlogy;
203 
204     const T z = (x - mu) / sigma;
205     if (xi != 0) {
206         return -log(sigma) + xlogy(-(cast(T) 1 / xi + 1), 1 + xi * z);
207     } else {
208         return -z;
209     }
210 }
211 
212 ///
213 version(mir_stat_test)
214 @safe pure nothrow @nogc
215 unittest {
216     import mir.math.common: log;
217     import mir.test: shouldApprox;
218 
219     1.0.generalizedParetoLPDF(1, 1, 0.5).shouldApprox == log(generalizedParetoPDF(1.0, 1, 1, 0.5));
220     2.0.generalizedParetoLPDF(1, 1, 0.5).shouldApprox == log(generalizedParetoPDF(2.0, 1, 1, 0.5));
221     3.0.generalizedParetoLPDF(2, 3, 0.25).shouldApprox == log(generalizedParetoPDF(3.0, 2, 3, 0.25));
222     5.0.generalizedParetoLPDF(2, 3, 0).shouldApprox == log(generalizedParetoPDF(5.0, 2, 3, 0));
223 }
The OpenD Programming Language
