tablepress / libraries / vendor / PhpSpreadsheet / Calculation / Statistical / Distributions / StudentT.php

StudentT.php in TablePress – Tables in WordPress made easy 3.4, at libraries/vendor/PhpSpreadsheet/Calculation/Statistical/Distributions/StudentT.php

333 lines 9.7 KB

Raw Download Zip

1	<?php
2
3	namespace TablePress\PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions;
4
5	use TablePress\PhpOffice\PhpSpreadsheet\Calculation\ArrayEnabled;
6	use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Exception;
7	use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError;
8
9	/**
10	* Some of this code is drived from Perl CPAN Statistical::Distributions.
11	* Its copyright statement is:
12	* Copyright 2003 Michael Kospach. All rights reserved.
13	*
14	* This library is free software; you can redistribute it and/or modify it under the same terms as Perl itself.
15	*/
16	class StudentT
17	{
18	use ArrayEnabled;
19
20	/**
21	* TDIST.
22	*
23	* Returns the probability of Student's T distribution.
24	*
25	* @param mixed $value Float value for the distribution
26	* Or can be an array of values
27	* @param mixed $degrees Integer value for degrees of freedom
28	* Or can be an array of values
29	* @param mixed $tails Integer value for the number of tails (1 or 2)
30	* Or can be an array of values
31	*
32	* @return array<mixed>\|float\|string The result, or a string containing an error
33	* If an array of numbers is passed as an argument, then the returned result will also be an array
34	* with the same dimensions
35	*/
36	public static function distribution($value, $degrees, $tails)
37	{
38	return self::calcDistribution($value, $degrees, $tails, \Closure::fromCallable([self::class, 'distribution']));
39	}
40
41	/**
42	* T.DIST.2T.
43	* Returns the two-tailed Student's t distribution.
44	*
45	* @return array<mixed>\|float\|string The result, or a string containing an error
46	* @param mixed $value
47	* @param mixed $degrees
48	*/
49	public static function tDotDistDot2T($value, $degrees)
50	{
51	return self::calcDistribution($value, $degrees, 2, \Closure::fromCallable([self::class, 'distribution']));
52	}
53
54	/**
55	* T.DIST.RT.
56	* Returns the right-tailed Student's t distribution.
57	*
58	* @return array<mixed>\|float\|string The result, or a string containing an error
59	* @param mixed $value
60	* @param mixed $degrees
61	*/
62	public static function tDotDistDotRT($value, $degrees)
63	{
64	return self::calcDistribution($value, $degrees, 1, \Closure::fromCallable([self::class, 'distribution']));
65	}
66
67	/**
68	* @return array<mixed>\|float\|string The result, or a string containing an error
69	* @param mixed $value
70	* @param mixed $degrees
71	* @param mixed $tails
72	*/
73	private static function calcDistribution($value, $degrees, $tails, callable $callback)
74	{
75	if (is_array($value) \|\| is_array($degrees) \|\| is_array($tails)) {
76	return self::evaluateArrayArguments($callback, $value, $degrees, $tails);
77	}
78
79	try {
80	$value = DistributionValidations::validateFloat($value);
81	$degrees = DistributionValidations::validateInt($degrees);
82	$tails = DistributionValidations::validateInt($tails);
83	} catch (Exception $e) {
84	return $e->getMessage();
85	}
86
87	if (($value < 0) \|\| ($degrees < 1) \|\| ($tails < 1) \|\| ($tails > 2)) {
88	return ExcelError::NAN();
89	}
90
91	return self::subTProb($value, $degrees, $tails);
92	}
93
94	/**
95	* Based on code from Perl CPAN Statistical::Distributions.
96	*/
97	private static function subTProb(float $x, int $n, int $tails): float
98	{
99	$w = atan2($x / sqrt($n), 1);
100	$z = cos($w) ** 2;
101	$y = 1;
102
103	for ($i = $n - 2; $i >= 2; $i -= 2) {
104	$y = 1 + ($i - 1) / $i * $z * $y;
105	}
106
107	if ($n % 2 == 0) {
108	$a = sin($w) / 2;
109	$b = 0.5;
110	} else {
111	$a = ($n == 1) ? 0 : (sin($w) * cos($w) / M_PI);
112	$b = 0.5 + $w / M_PI;
113	}
114
115	return $tails * max(0, 1 - $b - $a * $y);
116	}
117
118	/**
119	* T.DIST.
120	* Returns the Student's left-tailed t distribution,
121	* either as a cumulative distribution function (cdf) (TRUE)
122	* or as a probability density function (pdf) (FALSE),
123	* where TRUE/FALSE are the value of $cumulative parameter.
124	*
125	* "True" algoritm adapted from java.
126	* org.apache.commons.math3.distribution.TDistribution.
127	* "False" algorithm comes from:
128	* https://statproofbook.github.io/P/t-pdf.html
129	*
130	* @param mixed $cumulative Expecting bool. See above for explanation.
131	*
132	* @return array<mixed>\|float\|string The result, or a string containing an error
133	* @param mixed $value
134	* @param mixed $degrees
135	*/
136	public static function tDotDist($value, $degrees, $cumulative)
137	{
138	if (is_array($value) \|\| is_array($degrees) \|\| is_array($cumulative)) {
139	return self::evaluateArrayArguments(\Closure::fromCallable([self::class, 'tDotDist']), $value, $degrees, $cumulative);
140	}
141
142	try {
143	$value = DistributionValidations::validateFloat($value);
144	$degrees = DistributionValidations::validateInt($degrees);
145	$cumulative = DistributionValidations::validateBool($cumulative);
146	} catch (Exception $e) {
147	return $e->getMessage();
148	}
149
150	/** @var int $degrees */
151	if (($degrees < 1)) {
152	return ExcelError::NAN();
153	}
154	/** @var float $value */
155	if (!$cumulative) {
156	return self::tDotDistFalse($value, $degrees);
157	}
158
159	$f16 = $degrees / ($degrees + $value * $value);
160	$g16 = 0.5 * $degrees;
161	$h16 = 0.5;
162	$result = Beta::distribution($f16, $g16, $h16);
163	if (is_numeric($result)) {
164	$result = ($value < 0) ? (0.5 * $result) : (1 - 0.5 * $result);
165	}
166
167	return $result;
168	}
169
170	/**
171	* @return float\|string
172	*/
173	private static function tDotDistFalse(float $value, int $degrees)
174	{
175	$result = $k15 = Gamma::gamma(($degrees + 1) / 2);
176	if (is_numeric($k15)) {
177	$result = $k16 = Gamma::gamma($degrees / 2);
178	if (is_numeric($k16)) {
179	$k17 = sqrt(M_PI * $degrees);
180	$k18 = $k15 / ($k16 * $k17);
181	$k19 = $value * $value / $degrees + 1;
182	$k20 = -($degrees + 1) / 2;
183	$k21 = $k19 ** $k20;
184	$result = $k18 * $k21;
185	}
186	}
187
188	/** @var float\|string $result */
189	return $result;
190	}
191
192	/**
193	* TINV and T.INV.2T.
194	* Returns the two-tailed inverse of the Student t distribution.
195	*
196	* @param mixed $probability Float probability for the function
197	* Or can be an array of values
198	* @param mixed $degrees Integer value for degrees of freedom
199	* Or can be an array of values
200	*
201	* @return array<mixed>\|float\|string The result, or a string containing an error
202	* If an array of numbers is passed as an argument, then the returned result will also be an array
203	* with the same dimensions
204	*/
205	public static function inverse($probability, $degrees)
206	{
207	return self::calcInverse($probability, $degrees, 2, \Closure::fromCallable([self::class, 'inverse']));
208	}
209
210	/**
211	* @return array<mixed>\|float\|string The result, or a string containing an error
212	* @param mixed $probability
213	* @param mixed $degrees
214	*/
215	private static function calcInverse($probability, $degrees, int $tails, callable $callback2)
216	{
217	if (is_array($probability) \|\| is_array($degrees)) {
218	return self::evaluateArrayArguments($callback2, $probability, $degrees, $tails);
219	}
220
221	try {
222	$probability = DistributionValidations::validateProbability($probability);
223	$degrees = DistributionValidations::validateInt($degrees);
224	} catch (Exception $e) {
225	return $e->getMessage();
226	}
227
228	if ($degrees <= 0) {
229	return ExcelError::NAN();
230	}
231
232	$callback = fn ($value) => self::distribution($value, $degrees, $tails);
233
234	$newtonRaphson = new NewtonRaphson($callback);
235
236	$result = $newtonRaphson->execute($probability);
237	if (is_numeric($result) && $tails === 1) {
238	$result = -$result; // @codeCoverageIgnore
239	}
240
241	return $result;
242	}
243
244	/**
245	* T.INV.
246	* Returns the left-tailed inverse of the Student's t distribution.
247	*
248	* Based on code from Perl CPAN Statistical::Distributions.
249	*
250	* @return array<mixed>\|float\|string The result, or a string containing an error
251	* @param mixed $probability
252	* @param mixed $degrees
253	*/
254	public static function tDotInv($probability, $degrees)
255	{
256	if (is_array($probability) \|\| is_array($degrees)) {
257	return self::evaluateArrayArguments(\Closure::fromCallable([self::class, 'tDotInv']), $probability, $degrees);
258	}
259
260	try {
261	$probability = DistributionValidations::validateProbability($probability);
262	$degrees = DistributionValidations::validateInt($degrees);
263	} catch (Exception $e) {
264	return $e->getMessage();
265	}
266
267	if ($degrees < 1) {
268	return ExcelError::NAN();
269	}
270	if ($probability == 0.5) {
271	return 0.0;
272	}
273	if ($probability < 0.5) {
274	$result = self::tDotInv(1.0 - $probability, $degrees);
275
276	return is_numeric($result) ? -$result : $result;
277	}
278	$p = $probability;
279	$n = $degrees;
280	$u = self::subU($p);
281	$u2 = $u ** 2;
282
283	$a = ($u2 + 1) / 4;
284	$b = ((5 * $u2 + 16) * $u2 + 3) / 96;
285	$c = (((3 * $u2 + 19) * $u2 + 17) * $u2 - 15) / 384;
286	$d = ((((79 * $u2 + 776) * $u2 + 1482) * $u2 - 1920) * $u2 - 945) / 92160;
287	$e = (((((27 * $u2 + 339) * $u2 + 930) * $u2 - 1782) * $u2 - 765) * $u2 + 17955) / 368640;
288
289	$x = $u * (1 + ($a + ($b + ($c + ($d + $e / $n) / $n) / $n) / $n) / $n);
290
291	if ($n <= log10($p) ** 2 + 3) {
292	do {
293	$p1 = self::subTProb($x, $n, 1);
294	$n1 = $n + 1;
295	$delta = ($p1 - $p)
296	/ exp(($n1 * log($n1 / ($n + $x * $x))
297	+ log($n / $n1 / 2 / M_PI) - 1
298	+ (1 / $n1 - 1 / $n) / 6) / 2);
299	$x += $delta;
300	$round = sprintf('%.' . abs((int) (log10(abs($x)) - 4)) . 'F', $delta);
301	} while (($x) && ($round != 0));
302	}
303
304	return -$x;
305	}
306
307	/**
308	* Based on code from Perl CPAN Statistical::Distributions.
309	*/
310	private static function subU(float $p): float
311	{
312	$y = -log(4 * $p * (1 - $p));
313	$x = sqrt(
314	$y * (1.570796288
315	+ $y * (.03706987906
316	+ $y * (-.8364353589E-3
317	+ $y * (-.2250947176E-3
318	+ $y * (.6841218299E-5
319	+ $y * (0.5824238515E-5
320	+ $y * (-.104527497E-5
321	+ $y * (.8360937017E-7
322	+ $y * (-.3231081277E-8
323	+ $y * (.3657763036E-10
324	+ $y * .6936233982E-12))))))))))
325	);
326	if ($p > 0.5) {
327	$x = -$x;
328	}
329
330	return $x;
331	}
332	}
333