1<?php
2
3namespace PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions;
4
5use PhpOffice\PhpSpreadsheet\Calculation\Exception;
6use PhpOffice\PhpSpreadsheet\Calculation\Functions;
7
8class ChiSquared
9{
10    private const MAX_ITERATIONS = 256;
11
12    private const EPS = 2.22e-16;
13
14    /**
15     * CHIDIST.
16     *
17     * Returns the one-tailed probability of the chi-squared distribution.
18     *
19     * @param mixed $value Float value for which we want the probability
20     * @param mixed $degrees Integer degrees of freedom
21     *
22     * @return float|string
23     */
24    public static function distributionRightTail($value, $degrees)
25    {
26        $value = Functions::flattenSingleValue($value);
27        $degrees = Functions::flattenSingleValue($degrees);
28
29        try {
30            $value = DistributionValidations::validateFloat($value);
31            $degrees = DistributionValidations::validateInt($degrees);
32        } catch (Exception $e) {
33            return $e->getMessage();
34        }
35
36        if ($degrees < 1) {
37            return Functions::NAN();
38        }
39        if ($value < 0) {
40            if (Functions::getCompatibilityMode() == Functions::COMPATIBILITY_GNUMERIC) {
41                return 1;
42            }
43
44            return Functions::NAN();
45        }
46
47        return 1 - (Gamma::incompleteGamma($degrees / 2, $value / 2) / Gamma::gammaValue($degrees / 2));
48    }
49
50    /**
51     * CHIDIST.
52     *
53     * Returns the one-tailed probability of the chi-squared distribution.
54     *
55     * @param mixed $value Float value for which we want the probability
56     * @param mixed $degrees Integer degrees of freedom
57     * @param mixed $cumulative Boolean value indicating if we want the cdf (true) or the pdf (false)
58     *
59     * @return float|string
60     */
61    public static function distributionLeftTail($value, $degrees, $cumulative)
62    {
63        $value = Functions::flattenSingleValue($value);
64        $degrees = Functions::flattenSingleValue($degrees);
65        $cumulative = Functions::flattenSingleValue($cumulative);
66
67        try {
68            $value = DistributionValidations::validateFloat($value);
69            $degrees = DistributionValidations::validateInt($degrees);
70            $cumulative = DistributionValidations::validateBool($cumulative);
71        } catch (Exception $e) {
72            return $e->getMessage();
73        }
74
75        if ($degrees < 1) {
76            return Functions::NAN();
77        }
78        if ($value < 0) {
79            if (Functions::getCompatibilityMode() == Functions::COMPATIBILITY_GNUMERIC) {
80                return 1;
81            }
82
83            return Functions::NAN();
84        }
85
86        if ($cumulative === true) {
87            return 1 - self::distributionRightTail($value, $degrees);
88        }
89
90        return (($value ** (($degrees / 2) - 1) * exp(-$value / 2))) /
91            ((2 ** ($degrees / 2)) * Gamma::gammaValue($degrees / 2));
92    }
93
94    /**
95     * CHIINV.
96     *
97     * Returns the inverse of the right-tailed probability of the chi-squared distribution.
98     *
99     * @param mixed $probability Float probability at which you want to evaluate the distribution
100     * @param mixed $degrees Integer degrees of freedom
101     *
102     * @return float|string
103     */
104    public static function inverseRightTail($probability, $degrees)
105    {
106        $probability = Functions::flattenSingleValue($probability);
107        $degrees = Functions::flattenSingleValue($degrees);
108
109        try {
110            $probability = DistributionValidations::validateProbability($probability);
111            $degrees = DistributionValidations::validateInt($degrees);
112        } catch (Exception $e) {
113            return $e->getMessage();
114        }
115
116        if ($degrees < 1) {
117            return Functions::NAN();
118        }
119
120        $callback = function ($value) use ($degrees) {
121            return 1 - (Gamma::incompleteGamma($degrees / 2, $value / 2)
122                    / Gamma::gammaValue($degrees / 2));
123        };
124
125        $newtonRaphson = new NewtonRaphson($callback);
126
127        return $newtonRaphson->execute($probability);
128    }
129
130    /**
131     * CHIINV.
132     *
133     * Returns the inverse of the left-tailed probability of the chi-squared distribution.
134     *
135     * @param mixed $probability Float probability at which you want to evaluate the distribution
136     * @param mixed $degrees Integer degrees of freedom
137     *
138     * @return float|string
139     */
140    public static function inverseLeftTail($probability, $degrees)
141    {
142        $probability = Functions::flattenSingleValue($probability);
143        $degrees = Functions::flattenSingleValue($degrees);
144
145        try {
146            $probability = DistributionValidations::validateProbability($probability);
147            $degrees = DistributionValidations::validateInt($degrees);
148        } catch (Exception $e) {
149            return $e->getMessage();
150        }
151
152        if ($degrees < 1) {
153            return Functions::NAN();
154        }
155
156        return self::inverseLeftTailCalculation($probability, $degrees);
157    }
158
159    /**
160     * CHITEST.
161     *
162     * Uses the chi-square test to calculate the probability that the differences between two supplied data sets
163     *      (of observed and expected frequencies), are likely to be simply due to sampling error,
164     *      or if they are likely to be real.
165     *
166     * @param mixed $actual an array of observed frequencies
167     * @param mixed $expected an array of expected frequencies
168     *
169     * @return float|string
170     */
171    public static function test($actual, $expected)
172    {
173        $rows = count($actual);
174        $actual = Functions::flattenArray($actual);
175        $expected = Functions::flattenArray($expected);
176        $columns = count($actual) / $rows;
177
178        $countActuals = count($actual);
179        $countExpected = count($expected);
180        if ($countActuals !== $countExpected || $countActuals === 1) {
181            return Functions::NAN();
182        }
183
184        $result = 0.0;
185        for ($i = 0; $i < $countActuals; ++$i) {
186            if ($expected[$i] == 0.0) {
187                return Functions::DIV0();
188            } elseif ($expected[$i] < 0.0) {
189                return Functions::NAN();
190            }
191            $result += (($actual[$i] - $expected[$i]) ** 2) / $expected[$i];
192        }
193
194        $degrees = self::degrees($rows, $columns);
195
196        $result = self::distributionRightTail($result, $degrees);
197
198        return $result;
199    }
200
201    protected static function degrees(int $rows, int $columns): int
202    {
203        if ($rows === 1) {
204            return $columns - 1;
205        } elseif ($columns === 1) {
206            return $rows - 1;
207        }
208
209        return ($columns - 1) * ($rows - 1);
210    }
211
212    private static function inverseLeftTailCalculation(float $probability, int $degrees): float
213    {
214        // bracket the root
215        $min = 0;
216        $sd = sqrt(2.0 * $degrees);
217        $max = 2 * $sd;
218        $s = -1;
219
220        while ($s * self::pchisq($max, $degrees) > $probability * $s) {
221            $min = $max;
222            $max += 2 * $sd;
223        }
224
225        // Find root using bisection
226        $chi2 = 0.5 * ($min + $max);
227
228        while (($max - $min) > self::EPS * $chi2) {
229            if ($s * self::pchisq($chi2, $degrees) > $probability * $s) {
230                $min = $chi2;
231            } else {
232                $max = $chi2;
233            }
234            $chi2 = 0.5 * ($min + $max);
235        }
236
237        return $chi2;
238    }
239
240    private static function pchisq($chi2, $degrees)
241    {
242        return self::gammp($degrees, 0.5 * $chi2);
243    }
244
245    private static function gammp($n, $x)
246    {
247        if ($x < 0.5 * $n + 1) {
248            return self::gser($n, $x);
249        }
250
251        return 1 - self::gcf($n, $x);
252    }
253
254    // Return the incomplete gamma function P(n/2,x) evaluated by
255    // series representation. Algorithm from numerical recipe.
256    // Assume that n is a positive integer and x>0, won't check arguments.
257    // Relative error controlled by the eps parameter
258    private static function gser($n, $x)
259    {
260        $gln = Gamma::ln($n / 2);
261        $a = 0.5 * $n;
262        $ap = $a;
263        $sum = 1.0 / $a;
264        $del = $sum;
265        for ($i = 1; $i < 101; ++$i) {
266            ++$ap;
267            $del = $del * $x / $ap;
268            $sum += $del;
269            if ($del < $sum * self::EPS) {
270                break;
271            }
272        }
273
274        return $sum * exp(-$x + $a * log($x) - $gln);
275    }
276
277    // Return the incomplete gamma function Q(n/2,x) evaluated by
278    // its continued fraction representation. Algorithm from numerical recipe.
279    // Assume that n is a postive integer and x>0, won't check arguments.
280    // Relative error controlled by the eps parameter
281    private static function gcf($n, $x)
282    {
283        $gln = Gamma::ln($n / 2);
284        $a = 0.5 * $n;
285        $b = $x + 1 - $a;
286        $fpmin = 1.e-300;
287        $c = 1 / $fpmin;
288        $d = 1 / $b;
289        $h = $d;
290        for ($i = 1; $i < 101; ++$i) {
291            $an = -$i * ($i - $a);
292            $b += 2;
293            $d = $an * $d + $b;
294            if (abs($d) < $fpmin) {
295                $d = $fpmin;
296            }
297            $c = $b + $an / $c;
298            if (abs($c) < $fpmin) {
299                $c = $fpmin;
300            }
301            $d = 1 / $d;
302            $del = $d * $c;
303            $h = $h * $del;
304            if (abs($del - 1) < self::EPS) {
305                break;
306            }
307        }
308
309        return $h * exp(-$x + $a * log($x) - $gln);
310    }
311}
312