Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/stats/_ksstats.py: 9%
308 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1# Compute the two-sided one-sample Kolmogorov-Smirnov Prob(Dn <= d) where:
2# D_n = sup_x{|F_n(x) - F(x)|},
3# F_n(x) is the empirical CDF for a sample of size n {x_i: i=1,...,n},
4# F(x) is the CDF of a probability distribution.
5#
6# Exact methods:
7# Prob(D_n >= d) can be computed via a matrix algorithm of Durbin[1]
8# or a recursion algorithm due to Pomeranz[2].
9# Marsaglia, Tsang & Wang[3] gave a computation-efficient way to perform
10# the Durbin algorithm.
11# D_n >= d <==> D_n+ >= d or D_n- >= d (the one-sided K-S statistics), hence
12# Prob(D_n >= d) = 2*Prob(D_n+ >= d) - Prob(D_n+ >= d and D_n- >= d).
13# For d > 0.5, the latter intersection probability is 0.
14#
15# Approximate methods:
16# For d close to 0.5, ignoring that intersection term may still give a
17# reasonable approximation.
18# Li-Chien[4] and Korolyuk[5] gave an asymptotic formula extending
19# Kolmogorov's initial asymptotic, suitable for large d. (See
20# scipy.special.kolmogorov for that asymptotic)
21# Pelz-Good[6] used the functional equation for Jacobi theta functions to
22# transform the Li-Chien/Korolyuk formula produce a computational formula
23# suitable for small d.
24#
25# Simard and L'Ecuyer[7] provided an algorithm to decide when to use each of
26# the above approaches and it is that which is used here.
27#
28# Other approaches:
29# Carvalho[8] optimizes Durbin's matrix algorithm for large values of d.
30# Moscovich and Nadler[9] use FFTs to compute the convolutions.
32# References:
33# [1] Durbin J (1968).
34# "The Probability that the Sample Distribution Function Lies Between Two
35# Parallel Straight Lines."
36# Annals of Mathematical Statistics, 39, 398-411.
37# [2] Pomeranz J (1974).
38# "Exact Cumulative Distribution of the Kolmogorov-Smirnov Statistic for
39# Small Samples (Algorithm 487)."
40# Communications of the ACM, 17(12), 703-704.
41# [3] Marsaglia G, Tsang WW, Wang J (2003).
42# "Evaluating Kolmogorov's Distribution."
43# Journal of Statistical Software, 8(18), 1-4.
44# [4] LI-CHIEN, C. (1956).
45# "On the exact distribution of the statistics of A. N. Kolmogorov and
46# their asymptotic expansion."
47# Acta Matematica Sinica, 6, 55-81.
48# [5] KOROLYUK, V. S. (1960).
49# "Asymptotic analysis of the distribution of the maximum deviation in
50# the Bernoulli scheme."
51# Theor. Probability Appl., 4, 339-366.
52# [6] Pelz W, Good IJ (1976).
53# "Approximating the Lower Tail-areas of the Kolmogorov-Smirnov One-sample
54# Statistic."
55# Journal of the Royal Statistical Society, Series B, 38(2), 152-156.
56# [7] Simard, R., L'Ecuyer, P. (2011)
57# "Computing the Two-Sided Kolmogorov-Smirnov Distribution",
58# Journal of Statistical Software, Vol 39, 11, 1-18.
59# [8] Carvalho, Luis (2015)
60# "An Improved Evaluation of Kolmogorov's Distribution"
61# Journal of Statistical Software, Code Snippets; Vol 65(3), 1-8.
62# [9] Amit Moscovich, Boaz Nadler (2017)
63# "Fast calculation of boundary crossing probabilities for Poisson
64# processes",
65# Statistics & Probability Letters, Vol 123, 177-182.
68import numpy as np
69import scipy.special
70import scipy.special._ufuncs as scu
71from scipy._lib._finite_differences import _derivative
73_E128 = 128
74_EP128 = np.ldexp(np.longdouble(1), _E128)
75_EM128 = np.ldexp(np.longdouble(1), -_E128)
77_SQRT2PI = np.sqrt(2 * np.pi)
78_LOG_2PI = np.log(2 * np.pi)
79_MIN_LOG = -708
80_SQRT3 = np.sqrt(3)
81_PI_SQUARED = np.pi ** 2
82_PI_FOUR = np.pi ** 4
83_PI_SIX = np.pi ** 6
85# [Lifted from _loggamma.pxd.] If B_m are the Bernoulli numbers,
86# then Stirling coeffs are B_{2j}/(2j)/(2j-1) for j=8,...1.
87_STIRLING_COEFFS = [-2.955065359477124183e-2, 6.4102564102564102564e-3,
88 -1.9175269175269175269e-3, 8.4175084175084175084e-4,
89 -5.952380952380952381e-4, 7.9365079365079365079e-4,
90 -2.7777777777777777778e-3, 8.3333333333333333333e-2]
92def _log_nfactorial_div_n_pow_n(n):
93 # Computes n! / n**n
94 # = (n-1)! / n**(n-1)
95 # Uses Stirling's approximation, but removes n*log(n) up-front to
96 # avoid subtractive cancellation.
97 # = log(n)/2 - n + log(sqrt(2pi)) + sum B_{2j}/(2j)/(2j-1)/n**(2j-1)
98 rn = 1.0/n
99 return np.log(n)/2 - n + _LOG_2PI/2 + rn * np.polyval(_STIRLING_COEFFS, rn/n)
102def _clip_prob(p):
103 """clips a probability to range 0<=p<=1."""
104 return np.clip(p, 0.0, 1.0)
107def _select_and_clip_prob(cdfprob, sfprob, cdf=True):
108 """Selects either the CDF or SF, and then clips to range 0<=p<=1."""
109 p = np.where(cdf, cdfprob, sfprob)
110 return _clip_prob(p)
113def _kolmogn_DMTW(n, d, cdf=True):
114 r"""Computes the Kolmogorov CDF: Pr(D_n <= d) using the MTW approach to
115 the Durbin matrix algorithm.
117 Durbin (1968); Marsaglia, Tsang, Wang (2003). [1], [3].
118 """
119 # Write d = (k-h)/n, where k is positive integer and 0 <= h < 1
120 # Generate initial matrix H of size m*m where m=(2k-1)
121 # Compute k-th row of (n!/n^n) * H^n, scaling intermediate results.
122 # Requires memory O(m^2) and computation O(m^2 log(n)).
123 # Most suitable for small m.
125 if d >= 1.0:
126 return _select_and_clip_prob(1.0, 0.0, cdf)
127 nd = n * d
128 if nd <= 0.5:
129 return _select_and_clip_prob(0.0, 1.0, cdf)
130 k = int(np.ceil(nd))
131 h = k - nd
132 m = 2 * k - 1
134 H = np.zeros([m, m])
136 # Initialize: v is first column (and last row) of H
137 # v[j] = (1-h^(j+1)/(j+1)! (except for v[-1])
138 # w[j] = 1/(j)!
139 # q = k-th row of H (actually i!/n^i*H^i)
140 intm = np.arange(1, m + 1)
141 v = 1.0 - h ** intm
142 w = np.empty(m)
143 fac = 1.0
144 for j in intm:
145 w[j - 1] = fac
146 fac /= j # This might underflow. Isn't a problem.
147 v[j - 1] *= fac
148 tt = max(2 * h - 1.0, 0)**m - 2*h**m
149 v[-1] = (1.0 + tt) * fac
151 for i in range(1, m):
152 H[i - 1:, i] = w[:m - i + 1]
153 H[:, 0] = v
154 H[-1, :] = np.flip(v, axis=0)
156 Hpwr = np.eye(np.shape(H)[0]) # Holds intermediate powers of H
157 nn = n
158 expnt = 0 # Scaling of Hpwr
159 Hexpnt = 0 # Scaling of H
160 while nn > 0:
161 if nn % 2:
162 Hpwr = np.matmul(Hpwr, H)
163 expnt += Hexpnt
164 H = np.matmul(H, H)
165 Hexpnt *= 2
166 # Scale as needed.
167 if np.abs(H[k - 1, k - 1]) > _EP128:
168 H /= _EP128
169 Hexpnt += _E128
170 nn = nn // 2
172 p = Hpwr[k - 1, k - 1]
174 # Multiply by n!/n^n
175 for i in range(1, n + 1):
176 p = i * p / n
177 if np.abs(p) < _EM128:
178 p *= _EP128
179 expnt -= _E128
181 # unscale
182 if expnt != 0:
183 p = np.ldexp(p, expnt)
185 return _select_and_clip_prob(p, 1.0-p, cdf)
188def _pomeranz_compute_j1j2(i, n, ll, ceilf, roundf):
189 """Compute the endpoints of the interval for row i."""
190 if i == 0:
191 j1, j2 = -ll - ceilf - 1, ll + ceilf - 1
192 else:
193 # i + 1 = 2*ip1div2 + ip1mod2
194 ip1div2, ip1mod2 = divmod(i + 1, 2)
195 if ip1mod2 == 0: # i is odd
196 if ip1div2 == n + 1:
197 j1, j2 = n - ll - ceilf - 1, n + ll + ceilf - 1
198 else:
199 j1, j2 = ip1div2 - 1 - ll - roundf - 1, ip1div2 + ll - 1 + ceilf - 1
200 else:
201 j1, j2 = ip1div2 - 1 - ll - 1, ip1div2 + ll + roundf - 1
203 return max(j1 + 2, 0), min(j2, n)
206def _kolmogn_Pomeranz(n, x, cdf=True):
207 r"""Computes Pr(D_n <= d) using the Pomeranz recursion algorithm.
209 Pomeranz (1974) [2]
210 """
212 # V is n*(2n+2) matrix.
213 # Each row is convolution of the previous row and probabilities from a
214 # Poisson distribution.
215 # Desired CDF probability is n! V[n-1, 2n+1] (final entry in final row).
216 # Only two rows are needed at any given stage:
217 # - Call them V0 and V1.
218 # - Swap each iteration
219 # Only a few (contiguous) entries in each row can be non-zero.
220 # - Keep track of start and end (j1 and j2 below)
221 # - V0s and V1s track the start in the two rows
222 # Scale intermediate results as needed.
223 # Only a few different Poisson distributions can occur
224 t = n * x
225 ll = int(np.floor(t))
226 f = 1.0 * (t - ll) # fractional part of t
227 g = min(f, 1.0 - f)
228 ceilf = (1 if f > 0 else 0)
229 roundf = (1 if f > 0.5 else 0)
230 npwrs = 2 * (ll + 1) # Maximum number of powers needed in convolutions
231 gpower = np.empty(npwrs) # gpower = (g/n)^m/m!
232 twogpower = np.empty(npwrs) # twogpower = (2g/n)^m/m!
233 onem2gpower = np.empty(npwrs) # onem2gpower = ((1-2g)/n)^m/m!
234 # gpower etc are *almost* Poisson probs, just missing normalizing factor.
236 gpower[0] = 1.0
237 twogpower[0] = 1.0
238 onem2gpower[0] = 1.0
239 expnt = 0
240 g_over_n, two_g_over_n, one_minus_two_g_over_n = g/n, 2*g/n, (1 - 2*g)/n
241 for m in range(1, npwrs):
242 gpower[m] = gpower[m - 1] * g_over_n / m
243 twogpower[m] = twogpower[m - 1] * two_g_over_n / m
244 onem2gpower[m] = onem2gpower[m - 1] * one_minus_two_g_over_n / m
246 V0 = np.zeros([npwrs])
247 V1 = np.zeros([npwrs])
248 V1[0] = 1 # first row
249 V0s, V1s = 0, 0 # start indices of the two rows
251 j1, j2 = _pomeranz_compute_j1j2(0, n, ll, ceilf, roundf)
252 for i in range(1, 2 * n + 2):
253 # Preserve j1, V1, V1s, V0s from last iteration
254 k1 = j1
255 V0, V1 = V1, V0
256 V0s, V1s = V1s, V0s
257 V1.fill(0.0)
258 j1, j2 = _pomeranz_compute_j1j2(i, n, ll, ceilf, roundf)
259 if i == 1 or i == 2 * n + 1:
260 pwrs = gpower
261 else:
262 pwrs = (twogpower if i % 2 else onem2gpower)
263 ln2 = j2 - k1 + 1
264 if ln2 > 0:
265 conv = np.convolve(V0[k1 - V0s:k1 - V0s + ln2], pwrs[:ln2])
266 conv_start = j1 - k1 # First index to use from conv
267 conv_len = j2 - j1 + 1 # Number of entries to use from conv
268 V1[:conv_len] = conv[conv_start:conv_start + conv_len]
269 # Scale to avoid underflow.
270 if 0 < np.max(V1) < _EM128:
271 V1 *= _EP128
272 expnt -= _E128
273 V1s = V0s + j1 - k1
275 # multiply by n!
276 ans = V1[n - V1s]
277 for m in range(1, n + 1):
278 if np.abs(ans) > _EP128:
279 ans *= _EM128
280 expnt += _E128
281 ans *= m
283 # Undo any intermediate scaling
284 if expnt != 0:
285 ans = np.ldexp(ans, expnt)
286 ans = _select_and_clip_prob(ans, 1.0 - ans, cdf)
287 return ans
290def _kolmogn_PelzGood(n, x, cdf=True):
291 """Computes the Pelz-Good approximation to Prob(Dn <= x) with 0<=x<=1.
293 Start with Li-Chien, Korolyuk approximation:
294 Prob(Dn <= x) ~ K0(z) + K1(z)/sqrt(n) + K2(z)/n + K3(z)/n**1.5
295 where z = x*sqrt(n).
296 Transform each K_(z) using Jacobi theta functions into a form suitable
297 for small z.
298 Pelz-Good (1976). [6]
299 """
300 if x <= 0.0:
301 return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
302 if x >= 1.0:
303 return _select_and_clip_prob(1.0, 0.0, cdf=cdf)
305 z = np.sqrt(n) * x
306 zsquared, zthree, zfour, zsix = z**2, z**3, z**4, z**6
308 qlog = -_PI_SQUARED / 8 / zsquared
309 if qlog < _MIN_LOG: # z ~ 0.041743441416853426
310 return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
312 q = np.exp(qlog)
314 # Coefficients of terms in the sums for K1, K2 and K3
315 k1a = -zsquared
316 k1b = _PI_SQUARED / 4
318 k2a = 6 * zsix + 2 * zfour
319 k2b = (2 * zfour - 5 * zsquared) * _PI_SQUARED / 4
320 k2c = _PI_FOUR * (1 - 2 * zsquared) / 16
322 k3d = _PI_SIX * (5 - 30 * zsquared) / 64
323 k3c = _PI_FOUR * (-60 * zsquared + 212 * zfour) / 16
324 k3b = _PI_SQUARED * (135 * zfour - 96 * zsix) / 4
325 k3a = -30 * zsix - 90 * z**8
327 K0to3 = np.zeros(4)
328 # Use a Horner scheme to evaluate sum c_i q^(i^2)
329 # Reduces to a sum over odd integers.
330 maxk = int(np.ceil(16 * z / np.pi))
331 for k in range(maxk, 0, -1):
332 m = 2 * k - 1
333 msquared, mfour, msix = m**2, m**4, m**6
334 qpower = np.power(q, 8 * k)
335 coeffs = np.array([1.0,
336 k1a + k1b*msquared,
337 k2a + k2b*msquared + k2c*mfour,
338 k3a + k3b*msquared + k3c*mfour + k3d*msix])
339 K0to3 *= qpower
340 K0to3 += coeffs
341 K0to3 *= q
342 K0to3 *= _SQRT2PI
343 # z**10 > 0 as z > 0.04
344 K0to3 /= np.array([z, 6 * zfour, 72 * z**7, 6480 * z**10])
346 # Now do the other sum over the other terms, all integers k
347 # K_2: (pi^2 k^2) q^(k^2),
348 # K_3: (3pi^2 k^2 z^2 - pi^4 k^4)*q^(k^2)
349 # Don't expect much subtractive cancellation so use direct calculation
350 q = np.exp(-_PI_SQUARED / 2 / zsquared)
351 ks = np.arange(maxk, 0, -1)
352 ksquared = ks ** 2
353 sqrt3z = _SQRT3 * z
354 kspi = np.pi * ks
355 qpwers = q ** ksquared
356 k2extra = np.sum(ksquared * qpwers)
357 k2extra *= _PI_SQUARED * _SQRT2PI/(-36 * zthree)
358 K0to3[2] += k2extra
359 k3extra = np.sum((sqrt3z + kspi) * (sqrt3z - kspi) * ksquared * qpwers)
360 k3extra *= _PI_SQUARED * _SQRT2PI/(216 * zsix)
361 K0to3[3] += k3extra
362 powers_of_n = np.power(n * 1.0, np.arange(len(K0to3)) / 2.0)
363 K0to3 /= powers_of_n
365 if not cdf:
366 K0to3 *= -1
367 K0to3[0] += 1
369 Ksum = sum(K0to3)
370 return Ksum
373def _kolmogn(n, x, cdf=True):
374 """Computes the CDF(or SF) for the two-sided Kolmogorov-Smirnov statistic.
376 x must be of type float, n of type integer.
378 Simard & L'Ecuyer (2011) [7].
379 """
380 if np.isnan(n):
381 return n # Keep the same type of nan
382 if int(n) != n or n <= 0:
383 return np.nan
384 if x >= 1.0:
385 return _select_and_clip_prob(1.0, 0.0, cdf=cdf)
386 if x <= 0.0:
387 return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
388 t = n * x
389 if t <= 1.0: # Ruben-Gambino: 1/2n <= x <= 1/n
390 if t <= 0.5:
391 return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
392 if n <= 140:
393 prob = np.prod(np.arange(1, n+1) * (1.0/n) * (2*t - 1))
394 else:
395 prob = np.exp(_log_nfactorial_div_n_pow_n(n) + n * np.log(2*t-1))
396 return _select_and_clip_prob(prob, 1.0 - prob, cdf=cdf)
397 if t >= n - 1: # Ruben-Gambino
398 prob = 2 * (1.0 - x)**n
399 return _select_and_clip_prob(1 - prob, prob, cdf=cdf)
400 if x >= 0.5: # Exact: 2 * smirnov
401 prob = 2 * scipy.special.smirnov(n, x)
402 return _select_and_clip_prob(1.0 - prob, prob, cdf=cdf)
404 nxsquared = t * x
405 if n <= 140:
406 if nxsquared <= 0.754693:
407 prob = _kolmogn_DMTW(n, x, cdf=True)
408 return _select_and_clip_prob(prob, 1.0 - prob, cdf=cdf)
409 if nxsquared <= 4:
410 prob = _kolmogn_Pomeranz(n, x, cdf=True)
411 return _select_and_clip_prob(prob, 1.0 - prob, cdf=cdf)
412 # Now use Miller approximation of 2*smirnov
413 prob = 2 * scipy.special.smirnov(n, x)
414 return _select_and_clip_prob(1.0 - prob, prob, cdf=cdf)
416 # Split CDF and SF as they have different cutoffs on nxsquared.
417 if not cdf:
418 if nxsquared >= 370.0:
419 return 0.0
420 if nxsquared >= 2.2:
421 prob = 2 * scipy.special.smirnov(n, x)
422 return _clip_prob(prob)
423 # Fall through and compute the SF as 1.0-CDF
424 if nxsquared >= 18.0:
425 cdfprob = 1.0
426 elif n <= 100000 and n * x**1.5 <= 1.4:
427 cdfprob = _kolmogn_DMTW(n, x, cdf=True)
428 else:
429 cdfprob = _kolmogn_PelzGood(n, x, cdf=True)
430 return _select_and_clip_prob(cdfprob, 1.0 - cdfprob, cdf=cdf)
433def _kolmogn_p(n, x):
434 """Computes the PDF for the two-sided Kolmogorov-Smirnov statistic.
436 x must be of type float, n of type integer.
437 """
438 if np.isnan(n):
439 return n # Keep the same type of nan
440 if int(n) != n or n <= 0:
441 return np.nan
442 if x >= 1.0 or x <= 0:
443 return 0
444 t = n * x
445 if t <= 1.0:
446 # Ruben-Gambino: n!/n^n * (2t-1)^n -> 2 n!/n^n * n^2 * (2t-1)^(n-1)
447 if t <= 0.5:
448 return 0.0
449 if n <= 140:
450 prd = np.prod(np.arange(1, n) * (1.0 / n) * (2 * t - 1))
451 else:
452 prd = np.exp(_log_nfactorial_div_n_pow_n(n) + (n-1) * np.log(2 * t - 1))
453 return prd * 2 * n**2
454 if t >= n - 1:
455 # Ruben-Gambino : 1-2(1-x)**n -> 2n*(1-x)**(n-1)
456 return 2 * (1.0 - x) ** (n-1) * n
457 if x >= 0.5:
458 return 2 * scipy.stats.ksone.pdf(x, n)
460 # Just take a small delta.
461 # Ideally x +/- delta would stay within [i/n, (i+1)/n] for some integer a.
462 # as the CDF is a piecewise degree n polynomial.
463 # It has knots at 1/n, 2/n, ... (n-1)/n
464 # and is not a C-infinity function at the knots
465 delta = x / 2.0**16
466 delta = min(delta, x - 1.0/n)
467 delta = min(delta, 0.5 - x)
469 def _kk(_x):
470 return kolmogn(n, _x)
472 return _derivative(_kk, x, dx=delta, order=5)
475def _kolmogni(n, p, q):
476 """Computes the PPF/ISF of kolmogn.
478 n of type integer, n>= 1
479 p is the CDF, q the SF, p+q=1
480 """
481 if np.isnan(n):
482 return n # Keep the same type of nan
483 if int(n) != n or n <= 0:
484 return np.nan
485 if p <= 0:
486 return 1.0/n
487 if q <= 0:
488 return 1.0
489 delta = np.exp((np.log(p) - scipy.special.loggamma(n+1))/n)
490 if delta <= 1.0/n:
491 return (delta + 1.0 / n) / 2
492 x = -np.expm1(np.log(q/2.0)/n)
493 if x >= 1 - 1.0/n:
494 return x
495 x1 = scu._kolmogci(p)/np.sqrt(n)
496 x1 = min(x1, 1.0 - 1.0/n)
497 _f = lambda x: _kolmogn(n, x) - p
498 return scipy.optimize.brentq(_f, 1.0/n, x1, xtol=1e-14)
501def kolmogn(n, x, cdf=True):
502 """Computes the CDF for the two-sided Kolmogorov-Smirnov distribution.
504 The two-sided Kolmogorov-Smirnov distribution has as its CDF Pr(D_n <= x),
505 for a sample of size n drawn from a distribution with CDF F(t), where
506 D_n &= sup_t |F_n(t) - F(t)|, and
507 F_n(t) is the Empirical Cumulative Distribution Function of the sample.
509 Parameters
510 ----------
511 n : integer, array_like
512 the number of samples
513 x : float, array_like
514 The K-S statistic, float between 0 and 1
515 cdf : bool, optional
516 whether to compute the CDF(default=true) or the SF.
518 Returns
519 -------
520 cdf : ndarray
521 CDF (or SF it cdf is False) at the specified locations.
523 The return value has shape the result of numpy broadcasting n and x.
524 """
525 it = np.nditer([n, x, cdf, None],
526 op_dtypes=[None, np.float64, np.bool_, np.float64])
527 for _n, _x, _cdf, z in it:
528 if np.isnan(_n):
529 z[...] = _n
530 continue
531 if int(_n) != _n:
532 raise ValueError(f'n is not integral: {_n}')
533 z[...] = _kolmogn(int(_n), _x, cdf=_cdf)
534 result = it.operands[-1]
535 return result
538def kolmognp(n, x):
539 """Computes the PDF for the two-sided Kolmogorov-Smirnov distribution.
541 Parameters
542 ----------
543 n : integer, array_like
544 the number of samples
545 x : float, array_like
546 The K-S statistic, float between 0 and 1
548 Returns
549 -------
550 pdf : ndarray
551 The PDF at the specified locations
553 The return value has shape the result of numpy broadcasting n and x.
554 """
555 it = np.nditer([n, x, None])
556 for _n, _x, z in it:
557 if np.isnan(_n):
558 z[...] = _n
559 continue
560 if int(_n) != _n:
561 raise ValueError(f'n is not integral: {_n}')
562 z[...] = _kolmogn_p(int(_n), _x)
563 result = it.operands[-1]
564 return result
567def kolmogni(n, q, cdf=True):
568 """Computes the PPF(or ISF) for the two-sided Kolmogorov-Smirnov distribution.
570 Parameters
571 ----------
572 n : integer, array_like
573 the number of samples
574 q : float, array_like
575 Probabilities, float between 0 and 1
576 cdf : bool, optional
577 whether to compute the PPF(default=true) or the ISF.
579 Returns
580 -------
581 ppf : ndarray
582 PPF (or ISF if cdf is False) at the specified locations
584 The return value has shape the result of numpy broadcasting n and x.
585 """
586 it = np.nditer([n, q, cdf, None])
587 for _n, _q, _cdf, z in it:
588 if np.isnan(_n):
589 z[...] = _n
590 continue
591 if int(_n) != _n:
592 raise ValueError(f'n is not integral: {_n}')
593 _pcdf, _psf = (_q, 1-_q) if _cdf else (1-_q, _q)
594 z[...] = _kolmogni(int(_n), _pcdf, _psf)
595 result = it.operands[-1]
596 return result