mir.stat.distribution.cornish

1 /++
2 This module contains algorithms for the $(LINK2 https://en.wikipedia.org/wiki/Cornish%E2%80%93Fisher_expansion, Cornish-Fisher Expansion).
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.cornish_fisher;
13 
14 import mir.internal.utility: isFloatingPoint;
15 
16 /++
17 Approximates the inverse CDF of a continuous distribution using the Cornish-Fisher expansion.
18 
19 It is generally recommended to only use the Cornish-Fisher expansion with
20 distributions that are similar to the normal distribution. Extreme values of
21 `skewness` or `excessKurtosis` can result in poorer approximations.
22 
23 Params:
24     p = quantile to calculate inverse CDF
25     mu = mean
26     std = standard deviation
27     skewness = skewness
28     excessKurtosis = excess kurtosis (kurtosis - 3)
29 
30 See_also:
31     $(LINK2 https://en.wikipedia.org/wiki/Cornish%E2%80%93Fisher_expansion, Cornish-Fisher Expansion)
32 +/
33 T cornishFisherInvCDF(T)(const T p, const T mu, const T std, const T skewness, const T excessKurtosis)
34     if (isFloatingPoint!T)
35 {
36     return mu + std * cornishFisherInvCDF(p, skewness, excessKurtosis);
37 }
38 
39 ///
40 version(mir_stat_test)
41 @safe pure @nogc nothrow
42 unittest {
43     import mir.test: shouldApprox;
44 
45     0.99.cornishFisherInvCDF(0, 1, 0.1, 1).shouldApprox == 2.629904;
46     0.99.cornishFisherInvCDF(0.1, 0.2, 0.1, 1).shouldApprox == 0.6259808;
47 }
48 
49 /++
50 Ditto, but assumes mu = 0 and std = 1
51 
52 Params:
53     p = quantile to calculate inverse CDF
54     skewness = skewness (default = 0)
55     excessKurtosis = excess kurtosis (kurtosis - 3) (default = 0)
56 +/
57 T cornishFisherInvCDF(T)(const T p, const T skewness = 0, const T excessKurtosis = 0)
58     if (isFloatingPoint!T)
59 {
60     import mir.stat.distribution.normal: normalInvCDF;
61 
62     T x = normalInvCDF(p);
63     return x.cornishFisherInvCDFImpl(skewness, excessKurtosis);
64 }
65 
66 ///
67 version(mir_stat_test)
68 @safe pure @nogc nothrow
69 unittest {
70     import mir.test: shouldApprox;
71 
72     0.5.cornishFisherInvCDF.shouldApprox == 0;
73     0.5.cornishFisherInvCDF(1).shouldApprox == -0.1666667;
74     0.5.cornishFisherInvCDF(-1, 0).shouldApprox == 0.1666667;
75 
76     0.9.cornishFisherInvCDF(0.1, 0).shouldApprox == 1.292868;
77     0.9.cornishFisherInvCDF(0.1, 1).shouldApprox ==  1.220374;
78     0.9.cornishFisherInvCDF(0.1, -1).shouldApprox == 1.365363;
79 
80     0.99.cornishFisherInvCDF(0.1, 0).shouldApprox == 2.396116;
81     0.99.cornishFisherInvCDF(0.1, 1).shouldApprox == 2.629904;
82     0.99.cornishFisherInvCDF(0.1, -1).shouldApprox == 2.162328;
83 
84     0.01.cornishFisherInvCDF(0.1, 0).shouldApprox == -2.249053  ;
85     0.01.cornishFisherInvCDF(0.1, 1).shouldApprox == -2.482841;
86     0.01.cornishFisherInvCDF(0.1, -1).shouldApprox == -2.015265;
87 }
88 
89 package(mir.stat)
90 T cornishFisherInvCDFImpl(T)(const T x, const T skewness, const T excessKurtosis)
91     if (isFloatingPoint!T)
92 {
93     import mir.stat.distribution.normal: normalInvCDF;
94 
95     T x2 = x * x;
96     T x3 = x2 * x;
97     return x + (x2 - 1) * skewness / 6 + (x3 - 3 * x) * excessKurtosis / 24 - (2 * x3 - 5 * x) * skewness * skewness / 36;
98 }
99 
100 //
101 version(mir_stat_test)
102 @safe pure @nogc nothrow
103 unittest {
104     import mir.test: shouldApprox;
105 
106     0.0.cornishFisherInvCDFImpl(0, 0).shouldApprox == 0;
107     0.0.cornishFisherInvCDFImpl(1, 0).shouldApprox == -0.1666667;
108     0.0.cornishFisherInvCDFImpl(-1, 0).shouldApprox == 0.1666667;
109 
110     1.281552.cornishFisherInvCDFImpl(0.1, 0).shouldApprox == 1.292868;
111     1.281552.cornishFisherInvCDFImpl(0.1, 1).shouldApprox ==  1.220374;
112     1.281552.cornishFisherInvCDFImpl(0.1, -1).shouldApprox == 1.365363;
113 
114     2.326348.cornishFisherInvCDFImpl(0.1, 0).shouldApprox == 2.396116;
115     2.326348.cornishFisherInvCDFImpl(0.1, 1).shouldApprox == 2.629904;
116     2.326348.cornishFisherInvCDFImpl(0.1, -1).shouldApprox == 2.162328;
117 
118     (-2.326348).cornishFisherInvCDFImpl(0.1, 0).shouldApprox == -2.249053  ;
119     (-2.326348).cornishFisherInvCDFImpl(0.1, 1).shouldApprox == -2.482841;
120     (-2.326348).cornishFisherInvCDFImpl(0.1, -1).shouldApprox == -2.015265;
121 }
The OpenD Programming Language
