Code Coverage
 
Lines
Functions and Methods
Classes and Traits
Total
100.00% covered (success)
100.00%
38 / 38
100.00% covered (success)
100.00%
12 / 12
 CRAP
100.00% covered (success)
100.00%
1 / 1
ChiSquaredDistribution
100.00% covered (success)
100.00%
38 / 38
100.00% covered (success)
100.00%
12 / 12
 22
100.00% covered (success)
100.00%
1 / 1
 testHypothesis
100.00% covered (success)
100.00%
21 / 21
100.00% covered (success)
100.00%
1 / 1
 8
 getDegreesOfFreedom
100.00% covered (success)
100.00%
3 / 3
100.00% covered (success)
100.00%
1 / 1
 2
 getPdf
100.00% covered (success)
100.00%
3 / 3
100.00% covered (success)
100.00%
1 / 1
 2
 getCdf
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
 getMode
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
 getMean
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
 getMedian
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
 getVariance
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
 getStandardDeviation
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
 getMgf
100.00% covered (success)
100.00%
3 / 3
100.00% covered (success)
100.00%
1 / 1
 2
 getSkewness
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
 getExKurtosis
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
1<?php
2/**
3 * Jingga
4 *
5 * PHP Version 8.1
6 *
7 * @package   phpOMS\Math\Stochastic\Distribution
8 * @copyright Dennis Eichhorn
9 * @license   OMS License 2.0
10 * @version   1.0.0
11 * @link      https://jingga.app
12 */
13declare(strict_types=1);
14
15namespace phpOMS\Math\Stochastic\Distribution;
16
17use phpOMS\Math\Functions\Gamma;
18
19/**
20 * Chi squared distribution.
21 *
22 * Test if two variables are similarly distributed (goodness fit, test dependency of variables)
23 *
24 * @package phpOMS\Math\Stochastic\Distribution
25 * @license OMS License 2.0
26 * @link    https://jingga.app
27 * @since   1.0.0
28 */
29final class ChiSquaredDistribution
30{
31    /**
32     * Chi square table.
33     *
34     * @var array<int, array<string, float>>
35     * @since 1.0.0
36     */
37    public const TABLE = [
38        1   => ['0.995' => 0.000, '0.99' => 0.000, '0.975' => 0.001, '0.95' => 0.004, '0.90' => 0.016, '0.10' => 2.706, '0.05' => 3.841, '0.025' => 5.024, '0.01' => 6.635, '0.005' => 7.879],
39        2   => ['0.995' => 0.010, '0.99' => 0.020, '0.975' => 0.051, '0.95' => 0.103, '0.90' => 0.211, '0.10' => 4.605, '0.05' => 5.991, '0.025' => 7.378, '0.01' => 9.210, '0.005' => 10.597],
40        3   => ['0.995' => 0.072, '0.99' => 0.115, '0.975' => 0.216, '0.95' => 0.352, '0.90' => 0.584, '0.10' => 6.251, '0.05' => 7.815, '0.025' => 9.348, '0.01' => 11.345, '0.005' => 12.838],
41        4   => ['0.995' => 0.207, '0.99' => 0.297, '0.975' => 0.484, '0.95' => 0.711, '0.90' => 1.064, '0.10' => 7.779, '0.05' => 9.488, '0.025' => 11.143, '0.01' => 13.277, '0.005' => 14.860],
42        5   => ['0.995' => 0.412, '0.99' => 0.554, '0.975' => 0.831, '0.95' => 1.145, '0.90' => 1.610, '0.10' => 9.236, '0.05' => 11.070, '0.025' => 12.833, '0.01' => 15.086, '0.005' => 16.750],
43        6   => ['0.995' => 0.676, '0.99' => 0.872, '0.975' => 1.237, '0.95' => 1.635, '0.90' => 2.204, '0.10' => 10.645, '0.05' => 12.592, '0.025' => 14.449, '0.01' => 16.812, '0.005' => 18.548],
44        7   => ['0.995' => 0.989, '0.99' => 1.239, '0.975' => 1.690, '0.95' => 2.167, '0.90' => 2.833, '0.10' => 12.017, '0.05' => 14.067, '0.025' => 16.013, '0.01' => 18.475, '0.005' => 20.278],
45        8   => ['0.995' => 1.344, '0.99' => 1.646, '0.975' => 2.180, '0.95' => 2.733, '0.90' => 3.490, '0.10' => 13.362, '0.05' => 15.507, '0.025' => 17.535, '0.01' => 20.090, '0.005' => 21.955],
46        9   => ['0.995' => 1.735, '0.99' => 2.088, '0.975' => 2.700, '0.95' => 3.325, '0.90' => 4.168, '0.10' => 14.684, '0.05' => 16.919, '0.025' => 19.023, '0.01' => 21.666, '0.005' => 23.589],
47        10  => ['0.995' => 2.156, '0.99' => 2.558, '0.975' => 3.247, '0.95' => 3.940, '0.90' => 4.865, '0.10' => 15.987, '0.05' => 18.307, '0.025' => 20.483, '0.01' => 23.209, '0.005' => 25.188],
48        11  => ['0.995' => 2.603, '0.99' => 3.053, '0.975' => 3.816, '0.95' => 4.575, '0.90' => 5.578, '0.10' => 17.275, '0.05' => 19.675, '0.025' => 21.920, '0.01' => 24.725, '0.005' => 26.757],
49        12  => ['0.995' => 3.074, '0.99' => 3.571, '0.975' => 4.404, '0.95' => 5.226, '0.90' => 6.304, '0.10' => 18.549, '0.05' => 21.026, '0.025' => 23.337, '0.01' => 26.217, '0.005' => 28.300],
50        13  => ['0.995' => 3.565, '0.99' => 4.107, '0.975' => 5.009, '0.95' => 5.892, '0.90' => 7.042, '0.10' => 19.812, '0.05' => 22.362, '0.025' => 24.736, '0.01' => 27.688, '0.005' => 29.819],
51        14  => ['0.995' => 4.075, '0.99' => 4.660, '0.975' => 5.629, '0.95' => 6.571, '0.90' => 7.790, '0.10' => 21.064, '0.05' => 23.685, '0.025' => 26.119, '0.01' => 29.141, '0.005' => 31.319],
52        15  => ['0.995' => 4.601, '0.99' => 5.229, '0.975' => 6.262, '0.95' => 7.261, '0.90' => 8.547, '0.10' => 22.307, '0.05' => 24.996, '0.025' => 27.488, '0.01' => 30.578, '0.005' => 32.801],
53        16  => ['0.995' => 5.142, '0.99' => 5.812, '0.975' => 6.908, '0.95' => 7.962, '0.90' => 9.312, '0.10' => 23.542, '0.05' => 26.296, '0.025' => 28.845, '0.01' => 32.000, '0.005' => 34.267],
54        17  => ['0.995' => 5.697, '0.99' => 6.408, '0.975' => 7.564, '0.95' => 8.672, '0.90' => 10.085, '0.10' => 24.769, '0.05' => 27.587, '0.025' => 30.191, '0.01' => 33.409, '0.005' => 35.718],
55        18  => ['0.995' => 6.265, '0.99' => 7.015, '0.975' => 8.231, '0.95' => 9.390, '0.90' => 10.865, '0.10' => 25.989, '0.05' => 28.869, '0.025' => 31.526, '0.01' => 34.805, '0.005' => 37.156],
56        19  => ['0.995' => 6.844, '0.99' => 7.633, '0.975' => 8.907, '0.95' => 10.117, '0.90' => 11.651, '0.10' => 27.204, '0.05' => 30.144, '0.025' => 32.852, '0.01' => 36.191, '0.005' => 38.582],
57        20  => ['0.995' => 7.434, '0.99' => 8.260, '0.975' => 9.591, '0.95' => 10.851, '0.90' => 12.443, '0.10' => 28.412, '0.05' => 31.410, '0.025' => 34.170, '0.01' => 37.566, '0.005' => 39.997],
58        21  => ['0.995' => 8.034, '0.99' => 8.897, '0.975' => 10.283, '0.95' => 11.591, '0.90' => 13.240, '0.10' => 29.615, '0.05' => 32.671, '0.025' => 35.479, '0.01' => 38.932, '0.005' => 41.401],
59        22  => ['0.995' => 8.643, '0.99' => 9.542, '0.975' => 10.982, '0.95' => 12.338, '0.90' => 14.041, '0.10' => 30.813, '0.05' => 33.924, '0.025' => 36.781, '0.01' => 40.289, '0.005' => 42.796],
60        23  => ['0.995' => 9.260, '0.99' => 10.196, '0.975' => 11.689, '0.95' => 13.091, '0.90' => 14.848, '0.10' => 32.007, '0.05' => 35.172, '0.025' => 38.076, '0.01' => 41.638, '0.005' => 44.181],
61        24  => ['0.995' => 9.886, '0.99' => 10.856, '0.975' => 12.401, '0.95' => 13.848, '0.90' => 15.659, '0.10' => 33.196, '0.05' => 36.415, '0.025' => 39.364, '0.01' => 42.980, '0.005' => 45.559],
62        25  => ['0.995' => 10.520, '0.99' => 11.524, '0.975' => 13.120, '0.95' => 14.611, '0.90' => 16.473, '0.10' => 34.382, '0.05' => 37.652, '0.025' => 40.646, '0.01' => 44.314, '0.005' => 46.928],
63        26  => ['0.995' => 11.160, '0.99' => 12.198, '0.975' => 13.844, '0.95' => 15.379, '0.90' => 17.292, '0.10' => 35.563, '0.05' => 38.885, '0.025' => 41.923, '0.01' => 45.642, '0.005' => 48.290],
64        27  => ['0.995' => 11.808, '0.99' => 12.879, '0.975' => 14.573, '0.95' => 16.151, '0.90' => 18.114, '0.10' => 36.741, '0.05' => 40.113, '0.025' => 43.195, '0.01' => 46.963, '0.005' => 49.645],
65        28  => ['0.995' => 12.461, '0.99' => 13.565, '0.975' => 15.308, '0.95' => 16.928, '0.90' => 18.939, '0.10' => 37.916, '0.05' => 41.337, '0.025' => 44.461, '0.01' => 48.278, '0.005' => 50.993],
66        29  => ['0.995' => 13.121, '0.99' => 14.256, '0.975' => 16.047, '0.95' => 17.708, '0.90' => 19.768, '0.10' => 39.087, '0.05' => 42.557, '0.025' => 45.722, '0.01' => 49.588, '0.005' => 52.336],
67        30  => ['0.995' => 13.787, '0.99' => 14.953, '0.975' => 16.791, '0.95' => 18.493, '0.90' => 20.599, '0.10' => 40.256, '0.05' => 43.773, '0.025' => 46.979, '0.01' => 50.892, '0.005' => 53.672],
68        40  => ['0.995' => 20.707, '0.99' => 22.164, '0.975' => 24.433, '0.95' => 26.509, '0.90' => 29.051, '0.10' => 51.805, '0.05' => 55.758, '0.025' => 59.342, '0.01' => 63.691, '0.005' => 66.766],
69        50  => ['0.995' => 27.991, '0.99' => 29.707, '0.975' => 32.357, '0.95' => 34.764, '0.90' => 37.689, '0.10' => 63.167, '0.05' => 67.505, '0.025' => 71.420, '0.01' => 76.154, '0.005' => 79.490],
70        60  => ['0.995' => 35.534, '0.99' => 37.485, '0.975' => 40.482, '0.95' => 43.188, '0.90' => 46.459, '0.10' => 74.397, '0.05' => 79.082, '0.025' => 83.298, '0.01' => 88.379, '0.005' => 91.952],
71        70  => ['0.995' => 43.275, '0.99' => 45.442, '0.975' => 48.758, '0.95' => 51.739, '0.90' => 55.329, '0.10' => 85.527, '0.05' => 90.531, '0.025' => 95.023, '0.01' => 100.425, '0.005' => 104.215],
72        80  => ['0.995' => 51.172, '0.99' => 53.540, '0.975' => 57.153, '0.95' => 60.391, '0.90' => 64.278, '0.10' => 96.578, '0.05' => 101.879, '0.025' => 106.629, '0.01' => 112.329, '0.005' => 116.321],
73        90  => ['0.995' => 59.196, '0.99' => 61.754, '0.975' => 65.647, '0.95' => 69.126, '0.90' => 73.291, '0.10' => 107.565, '0.05' => 113.145, '0.025' => 118.136, '0.01' => 124.116, '0.005' => 128.299],
74        100 => ['0.995' => 67.328, '0.99' => 70.065, '0.975' => 74.222, '0.95' => 77.929, '0.90' => 82.358, '0.10' => 118.498, '0.05' => 124.342, '0.025' => 129.561, '0.01' => 135.807, '0.005' => 140.169],
75    ];
76
77    /**
78     * Test hypthesis.
79     *
80     * Goodness of fit test.
81     *
82     * @param array $dataset      Values
83     * @param array $expected     Expected values based on probability
84     * @param float $significance Significance
85     * @param int   $df           Degrees of freedom (optional)
86     *
87     * @return array
88     *
89     * @throws \Exception
90     *
91     * @since 1.0.0
92     */
93    public static function testHypothesis(array $dataset, array $expected, float $significance = 0.05, int $df = 0) : array
94    {
95        if (($count = \count($dataset)) !== \count($expected)) {
96            throw new \Exception('Dimension');
97        }
98
99        $sum = 0.0;
100
101        for ($i = 0; $i < $count; ++$i) {
102            $sum += ($dataset[$i] - $expected[$i]) ** 2 / $expected[$i];
103        }
104
105        $p = null;
106
107        if ($df === 0) {
108            $df = self::getDegreesOfFreedom($dataset);
109        }
110
111        if (!\defined('self::TABLE') || !\array_key_exists($df, self::TABLE)) {
112            throw new \Exception('Degrees of freedom not supported');
113        }
114
115        foreach (self::TABLE[$df] as $key => $value) {
116            if ($value > $sum) {
117                $p = $key;
118                break;
119            }
120        }
121
122        $p = $p ?? 0;
123
124        return [
125            'Chi2' => $sum,
126            'P'    => $p,
127            'H0'   => ($p > $significance),
128            'df'   => $df,
129        ];
130    }
131
132    /**
133     * Get degrees of freedom of array.
134     *
135     * @param array $values Value matrix or vector (N or NxM)
136     *
137     * @return int
138     *
139     * @since 1.0.0
140     */
141    public static function getDegreesOfFreedom(array $values) : int
142    {
143        if (\is_array($first = \reset($values))) {
144            return (\count($values) - 1) * (\count($first) - 1);
145        } else {
146            return \count($values) - 1;
147        }
148    }
149
150    /**
151     * Get probability density function.
152     *
153     * @param float $x  Value x
154     * @param int   $df Degreegs of freedom
155     *
156     * @return float
157     *
158     * @throws \OutOfBoundsException
159     *
160     * @since 1.0.0
161     */
162    public static function getPdf(float $x, int $df) : float
163    {
164        if ($x < 0) {
165            throw new \OutOfBoundsException('Out of bounds');
166        }
167
168        return 1 / (\pow(2, $df / 2) * Gamma::gamma($df / 2)) * \pow($x, $df / 2 - 1) * \exp(-$x / 2);
169    }
170
171    /**
172     * Get cumulative density function.
173     *
174     * @param float $x  Value x
175     * @param int   $df Degreegs of freedom
176     *
177     * @return float
178     *
179     * @since 1.0.0
180     */
181    public static function getCdf(float $x, int $df) : float
182    {
183        return 1 / Gamma::gamma($df / 2) * Gamma::incompleteGammaFirst($df / 2, $x / 2);
184    }
185
186    /**
187     * Get mode.
188     *
189     * @param int $df Degrees of freedom
190     *
191     * @return int
192     *
193     * @since 1.0.0
194     */
195    public static function getMode(int $df) : int
196    {
197        return \max($df - 2, 0);
198    }
199
200    /**
201     * Get expected value.
202     *
203     * @param int $df Degrees of freedom
204     *
205     * @return float
206     *
207     * @since 1.0.0
208     */
209    public static function getMean(int $df) : float
210    {
211        return $df;
212    }
213
214    /**
215     * Get median.
216     *
217     * @param int $df Degrees of freedom
218     *
219     * @return float
220     *
221     * @since 1.0.0
222     */
223    public static function getMedian(int $df) : float
224    {
225        return $df * (1 - 2 / (9 * $df)) ** 3;
226    }
227
228    /**
229     * Get variance.
230     *
231     * @param int $df Degrees of freedom
232     *
233     * @return float
234     *
235     * @since 1.0.0
236     */
237    public static function getVariance(int $df) : float
238    {
239        return 2 * $df;
240    }
241
242    /**
243     * Get standard deviation.
244     *
245     * @param int $df Degrees of freedom
246     *
247     * @return float
248     *
249     * @since 1.0.0
250     */
251    public static function getStandardDeviation(int $df) : float
252    {
253        return \sqrt(self::getVariance($df));
254    }
255
256    /**
257     * Get moment generating function.
258     *
259     * @param int   $df Degrees of freedom
260     * @param float $t  Value t
261     *
262     * @return float
263     *
264     * @throws \OutOfBoundsException
265     *
266     * @since 1.0.0
267     */
268    public static function getMgf(int $df, float $t) : float
269    {
270        if ($t > 0.5) {
271            throw new \OutOfBoundsException('Out of bounds');
272        }
273
274        return \pow(1 - 2 * $t, -$df / 2);
275    }
276
277    /**
278     * Get skewness.
279     *
280     * @param int $df Degrees of freedom
281     *
282     * @return float
283     *
284     * @since 1.0.0
285     */
286    public static function getSkewness(int $df) : float
287    {
288        return \sqrt(8 / $df);
289    }
290
291    /**
292     * Get Ex. kurtosis.
293     *
294     * @param int $df Degrees of freedom
295     *
296     * @return float
297     *
298     * @since 1.0.0
299     */
300    public static function getExKurtosis(int $df) : float
301    {
302        return 12 / $df;
303    }
304}