Coverage for /usr/lib/python3/dist-packages/scipy/stats/_ksstats.py: 9%
309 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
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]
93def _log_nfactorial_div_n_pow_n(n):
94 # Computes n! / n**n
95 # = (n-1)! / n**(n-1)
96 # Uses Stirling's approximation, but removes n*log(n) up-front to
97 # avoid subtractive cancellation.
98 # = log(n)/2 - n + log(sqrt(2pi)) + sum B_{2j}/(2j)/(2j-1)/n**(2j-1)
99 rn = 1.0/n
100 return np.log(n)/2 - n + _LOG_2PI/2 + rn * np.polyval(_STIRLING_COEFFS, rn/n)
103def _clip_prob(p):
104 """clips a probability to range 0<=p<=1."""
105 return np.clip(p, 0.0, 1.0)
108def _select_and_clip_prob(cdfprob, sfprob, cdf=True):
109 """Selects either the CDF or SF, and then clips to range 0<=p<=1."""
110 p = np.where(cdf, cdfprob, sfprob)
111 return _clip_prob(p)
114def _kolmogn_DMTW(n, d, cdf=True):
115 r"""Computes the Kolmogorov CDF: Pr(D_n <= d) using the MTW approach to
116 the Durbin matrix algorithm.
118 Durbin (1968); Marsaglia, Tsang, Wang (2003). [1], [3].
119 """
120 # Write d = (k-h)/n, where k is positive integer and 0 <= h < 1
121 # Generate initial matrix H of size m*m where m=(2k-1)
122 # Compute k-th row of (n!/n^n) * H^n, scaling intermediate results.
123 # Requires memory O(m^2) and computation O(m^2 log(n)).
124 # Most suitable for small m.
126 if d >= 1.0:
127 return _select_and_clip_prob(1.0, 0.0, cdf)
128 nd = n * d
129 if nd <= 0.5:
130 return _select_and_clip_prob(0.0, 1.0, cdf)
131 k = int(np.ceil(nd))
132 h = k - nd
133 m = 2 * k - 1
135 H = np.zeros([m, m])
137 # Initialize: v is first column (and last row) of H
138 # v[j] = (1-h^(j+1)/(j+1)! (except for v[-1])
139 # w[j] = 1/(j)!
140 # q = k-th row of H (actually i!/n^i*H^i)
141 intm = np.arange(1, m + 1)
142 v = 1.0 - h ** intm
143 w = np.empty(m)
144 fac = 1.0
145 for j in intm:
146 w[j - 1] = fac
147 fac /= j # This might underflow. Isn't a problem.
148 v[j - 1] *= fac
149 tt = max(2 * h - 1.0, 0)**m - 2*h**m
150 v[-1] = (1.0 + tt) * fac
152 for i in range(1, m):
153 H[i - 1:, i] = w[:m - i + 1]
154 H[:, 0] = v
155 H[-1, :] = np.flip(v, axis=0)
157 Hpwr = np.eye(np.shape(H)[0]) # Holds intermediate powers of H
158 nn = n
159 expnt = 0 # Scaling of Hpwr
160 Hexpnt = 0 # Scaling of H
161 while nn > 0:
162 if nn % 2:
163 Hpwr = np.matmul(Hpwr, H)
164 expnt += Hexpnt
165 H = np.matmul(H, H)
166 Hexpnt *= 2
167 # Scale as needed.
168 if np.abs(H[k - 1, k - 1]) > _EP128:
169 H /= _EP128
170 Hexpnt += _E128
171 nn = nn // 2
173 p = Hpwr[k - 1, k - 1]
175 # Multiply by n!/n^n
176 for i in range(1, n + 1):
177 p = i * p / n
178 if np.abs(p) < _EM128:
179 p *= _EP128
180 expnt -= _E128
182 # unscale
183 if expnt != 0:
184 p = np.ldexp(p, expnt)
186 return _select_and_clip_prob(p, 1.0-p, cdf)
189def _pomeranz_compute_j1j2(i, n, ll, ceilf, roundf):
190 """Compute the endpoints of the interval for row i."""
191 if i == 0:
192 j1, j2 = -ll - ceilf - 1, ll + ceilf - 1
193 else:
194 # i + 1 = 2*ip1div2 + ip1mod2
195 ip1div2, ip1mod2 = divmod(i + 1, 2)
196 if ip1mod2 == 0: # i is odd
197 if ip1div2 == n + 1:
198 j1, j2 = n - ll - ceilf - 1, n + ll + ceilf - 1
199 else:
200 j1, j2 = ip1div2 - 1 - ll - roundf - 1, ip1div2 + ll - 1 + ceilf - 1
201 else:
202 j1, j2 = ip1div2 - 1 - ll - 1, ip1div2 + ll + roundf - 1
204 return max(j1 + 2, 0), min(j2, n)
207def _kolmogn_Pomeranz(n, x, cdf=True):
208 r"""Computes Pr(D_n <= d) using the Pomeranz recursion algorithm.
210 Pomeranz (1974) [2]
211 """
213 # V is n*(2n+2) matrix.
214 # Each row is convolution of the previous row and probabilities from a
215 # Poisson distribution.
216 # Desired CDF probability is n! V[n-1, 2n+1] (final entry in final row).
217 # Only two rows are needed at any given stage:
218 # - Call them V0 and V1.
219 # - Swap each iteration
220 # Only a few (contiguous) entries in each row can be non-zero.
221 # - Keep track of start and end (j1 and j2 below)
222 # - V0s and V1s track the start in the two rows
223 # Scale intermediate results as needed.
224 # Only a few different Poisson distributions can occur
225 t = n * x
226 ll = int(np.floor(t))
227 f = 1.0 * (t - ll) # fractional part of t
228 g = min(f, 1.0 - f)
229 ceilf = (1 if f > 0 else 0)
230 roundf = (1 if f > 0.5 else 0)
231 npwrs = 2 * (ll + 1) # Maximum number of powers needed in convolutions
232 gpower = np.empty(npwrs) # gpower = (g/n)^m/m!
233 twogpower = np.empty(npwrs) # twogpower = (2g/n)^m/m!
234 onem2gpower = np.empty(npwrs) # onem2gpower = ((1-2g)/n)^m/m!
235 # gpower etc are *almost* Poisson probs, just missing normalizing factor.
237 gpower[0] = 1.0
238 twogpower[0] = 1.0
239 onem2gpower[0] = 1.0
240 expnt = 0
241 g_over_n, two_g_over_n, one_minus_two_g_over_n = g/n, 2*g/n, (1 - 2*g)/n
242 for m in range(1, npwrs):
243 gpower[m] = gpower[m - 1] * g_over_n / m
244 twogpower[m] = twogpower[m - 1] * two_g_over_n / m
245 onem2gpower[m] = onem2gpower[m - 1] * one_minus_two_g_over_n / m
247 V0 = np.zeros([npwrs])
248 V1 = np.zeros([npwrs])
249 V1[0] = 1 # first row
250 V0s, V1s = 0, 0 # start indices of the two rows
252 j1, j2 = _pomeranz_compute_j1j2(0, n, ll, ceilf, roundf)
253 for i in range(1, 2 * n + 2):
254 # Preserve j1, V1, V1s, V0s from last iteration
255 k1 = j1
256 V0, V1 = V1, V0
257 V0s, V1s = V1s, V0s
258 V1.fill(0.0)
259 j1, j2 = _pomeranz_compute_j1j2(i, n, ll, ceilf, roundf)
260 if i == 1 or i == 2 * n + 1:
261 pwrs = gpower
262 else:
263 pwrs = (twogpower if i % 2 else onem2gpower)
264 ln2 = j2 - k1 + 1
265 if ln2 > 0:
266 conv = np.convolve(V0[k1 - V0s:k1 - V0s + ln2], pwrs[:ln2])
267 conv_start = j1 - k1 # First index to use from conv
268 conv_len = j2 - j1 + 1 # Number of entries to use from conv
269 V1[:conv_len] = conv[conv_start:conv_start + conv_len]
270 # Scale to avoid underflow.
271 if 0 < np.max(V1) < _EM128:
272 V1 *= _EP128
273 expnt -= _E128
274 V1s = V0s + j1 - k1
276 # multiply by n!
277 ans = V1[n - V1s]
278 for m in range(1, n + 1):
279 if np.abs(ans) > _EP128:
280 ans *= _EM128
281 expnt += _E128
282 ans *= m
284 # Undo any intermediate scaling
285 if expnt != 0:
286 ans = np.ldexp(ans, expnt)
287 ans = _select_and_clip_prob(ans, 1.0 - ans, cdf)
288 return ans
291def _kolmogn_PelzGood(n, x, cdf=True):
292 """Computes the Pelz-Good approximation to Prob(Dn <= x) with 0<=x<=1.
294 Start with Li-Chien, Korolyuk approximation:
295 Prob(Dn <= x) ~ K0(z) + K1(z)/sqrt(n) + K2(z)/n + K3(z)/n**1.5
296 where z = x*sqrt(n).
297 Transform each K_(z) using Jacobi theta functions into a form suitable
298 for small z.
299 Pelz-Good (1976). [6]
300 """
301 if x <= 0.0:
302 return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
303 if x >= 1.0:
304 return _select_and_clip_prob(1.0, 0.0, cdf=cdf)
306 z = np.sqrt(n) * x
307 zsquared, zthree, zfour, zsix = z**2, z**3, z**4, z**6
309 qlog = -_PI_SQUARED / 8 / zsquared
310 if qlog < _MIN_LOG: # z ~ 0.041743441416853426
311 return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
313 q = np.exp(qlog)
315 # Coefficients of terms in the sums for K1, K2 and K3
316 k1a = -zsquared
317 k1b = _PI_SQUARED / 4
319 k2a = 6 * zsix + 2 * zfour
320 k2b = (2 * zfour - 5 * zsquared) * _PI_SQUARED / 4
321 k2c = _PI_FOUR * (1 - 2 * zsquared) / 16
323 k3d = _PI_SIX * (5 - 30 * zsquared) / 64
324 k3c = _PI_FOUR * (-60 * zsquared + 212 * zfour) / 16
325 k3b = _PI_SQUARED * (135 * zfour - 96 * zsix) / 4
326 k3a = -30 * zsix - 90 * z**8
328 K0to3 = np.zeros(4)
329 # Use a Horner scheme to evaluate sum c_i q^(i^2)
330 # Reduces to a sum over odd integers.
331 maxk = int(np.ceil(16 * z / np.pi))
332 for k in range(maxk, 0, -1):
333 m = 2 * k - 1
334 msquared, mfour, msix = m**2, m**4, m**6
335 qpower = np.power(q, 8 * k)
336 coeffs = np.array([1.0,
337 k1a + k1b*msquared,
338 k2a + k2b*msquared + k2c*mfour,
339 k3a + k3b*msquared + k3c*mfour + k3d*msix])
340 K0to3 *= qpower
341 K0to3 += coeffs
342 K0to3 *= q
343 K0to3 *= _SQRT2PI
344 # z**10 > 0 as z > 0.04
345 K0to3 /= np.array([z, 6 * zfour, 72 * z**7, 6480 * z**10])
347 # Now do the other sum over the other terms, all integers k
348 # K_2: (pi^2 k^2) q^(k^2),
349 # K_3: (3pi^2 k^2 z^2 - pi^4 k^4)*q^(k^2)
350 # Don't expect much subtractive cancellation so use direct calculation
351 q = np.exp(-_PI_SQUARED / 2 / zsquared)
352 ks = np.arange(maxk, 0, -1)
353 ksquared = ks ** 2
354 sqrt3z = _SQRT3 * z
355 kspi = np.pi * ks
356 qpwers = q ** ksquared
357 k2extra = np.sum(ksquared * qpwers)
358 k2extra *= _PI_SQUARED * _SQRT2PI/(-36 * zthree)
359 K0to3[2] += k2extra
360 k3extra = np.sum((sqrt3z + kspi) * (sqrt3z - kspi) * ksquared * qpwers)
361 k3extra *= _PI_SQUARED * _SQRT2PI/(216 * zsix)
362 K0to3[3] += k3extra
363 powers_of_n = np.power(n * 1.0, np.arange(len(K0to3)) / 2.0)
364 K0to3 /= powers_of_n
366 if not cdf:
367 K0to3 *= -1
368 K0to3[0] += 1
370 Ksum = sum(K0to3)
371 return Ksum
374def _kolmogn(n, x, cdf=True):
375 """Computes the CDF(or SF) for the two-sided Kolmogorov-Smirnov statistic.
377 x must be of type float, n of type integer.
379 Simard & L'Ecuyer (2011) [7].
380 """
381 if np.isnan(n):
382 return n # Keep the same type of nan
383 if int(n) != n or n <= 0:
384 return np.nan
385 if x >= 1.0:
386 return _select_and_clip_prob(1.0, 0.0, cdf=cdf)
387 if x <= 0.0:
388 return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
389 t = n * x
390 if t <= 1.0: # Ruben-Gambino: 1/2n <= x <= 1/n
391 if t <= 0.5:
392 return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
393 if n <= 140:
394 prob = np.prod(np.arange(1, n+1) * (1.0/n) * (2*t - 1))
395 else:
396 prob = np.exp(_log_nfactorial_div_n_pow_n(n) + n * np.log(2*t-1))
397 return _select_and_clip_prob(prob, 1.0 - prob, cdf=cdf)
398 if t >= n - 1: # Ruben-Gambino
399 prob = 2 * (1.0 - x)**n
400 return _select_and_clip_prob(1 - prob, prob, cdf=cdf)
401 if x >= 0.5: # Exact: 2 * smirnov
402 prob = 2 * scipy.special.smirnov(n, x)
403 return _select_and_clip_prob(1.0 - prob, prob, cdf=cdf)
405 nxsquared = t * x
406 if n <= 140:
407 if nxsquared <= 0.754693:
408 prob = _kolmogn_DMTW(n, x, cdf=True)
409 return _select_and_clip_prob(prob, 1.0 - prob, cdf=cdf)
410 if nxsquared <= 4:
411 prob = _kolmogn_Pomeranz(n, x, cdf=True)
412 return _select_and_clip_prob(prob, 1.0 - prob, cdf=cdf)
413 # Now use Miller approximation of 2*smirnov
414 prob = 2 * scipy.special.smirnov(n, x)
415 return _select_and_clip_prob(1.0 - prob, prob, cdf=cdf)
417 # Split CDF and SF as they have different cutoffs on nxsquared.
418 if not cdf:
419 if nxsquared >= 370.0:
420 return 0.0
421 if nxsquared >= 2.2:
422 prob = 2 * scipy.special.smirnov(n, x)
423 return _clip_prob(prob)
424 # Fall through and compute the SF as 1.0-CDF
425 if nxsquared >= 18.0:
426 cdfprob = 1.0
427 elif n <= 100000 and n * x**1.5 <= 1.4:
428 cdfprob = _kolmogn_DMTW(n, x, cdf=True)
429 else:
430 cdfprob = _kolmogn_PelzGood(n, x, cdf=True)
431 return _select_and_clip_prob(cdfprob, 1.0 - cdfprob, cdf=cdf)
434def _kolmogn_p(n, x):
435 """Computes the PDF for the two-sided Kolmogorov-Smirnov statistic.
437 x must be of type float, n of type integer.
438 """
439 if np.isnan(n):
440 return n # Keep the same type of nan
441 if int(n) != n or n <= 0:
442 return np.nan
443 if x >= 1.0 or x <= 0:
444 return 0
445 t = n * x
446 if t <= 1.0:
447 # Ruben-Gambino: n!/n^n * (2t-1)^n -> 2 n!/n^n * n^2 * (2t-1)^(n-1)
448 if t <= 0.5:
449 return 0.0
450 if n <= 140:
451 prd = np.prod(np.arange(1, n) * (1.0 / n) * (2 * t - 1))
452 else:
453 prd = np.exp(_log_nfactorial_div_n_pow_n(n) + (n-1) * np.log(2 * t - 1))
454 return prd * 2 * n**2
455 if t >= n - 1:
456 # Ruben-Gambino : 1-2(1-x)**n -> 2n*(1-x)**(n-1)
457 return 2 * (1.0 - x) ** (n-1) * n
458 if x >= 0.5:
459 return 2 * scipy.stats.ksone.pdf(x, n)
461 # Just take a small delta.
462 # Ideally x +/- delta would stay within [i/n, (i+1)/n] for some integer a.
463 # as the CDF is a piecewise degree n polynomial.
464 # It has knots at 1/n, 2/n, ... (n-1)/n
465 # and is not a C-infinity function at the knots
466 delta = x / 2.0**16
467 delta = min(delta, x - 1.0/n)
468 delta = min(delta, 0.5 - x)
470 def _kk(_x):
471 return kolmogn(n, _x)
473 return _derivative(_kk, x, dx=delta, order=5)
476def _kolmogni(n, p, q):
477 """Computes the PPF/ISF of kolmogn.
479 n of type integer, n>= 1
480 p is the CDF, q the SF, p+q=1
481 """
482 if np.isnan(n):
483 return n # Keep the same type of nan
484 if int(n) != n or n <= 0:
485 return np.nan
486 if p <= 0:
487 return 1.0/n
488 if q <= 0:
489 return 1.0
490 delta = np.exp((np.log(p) - scipy.special.loggamma(n+1))/n)
491 if delta <= 1.0/n:
492 return (delta + 1.0 / n) / 2
493 x = -np.expm1(np.log(q/2.0)/n)
494 if x >= 1 - 1.0/n:
495 return x
496 x1 = scu._kolmogci(p)/np.sqrt(n)
497 x1 = min(x1, 1.0 - 1.0/n)
499 def _f(x):
500 return _kolmogn(n, x) - p
502 return scipy.optimize.brentq(_f, 1.0/n, x1, xtol=1e-14)
505def kolmogn(n, x, cdf=True):
506 """Computes the CDF for the two-sided Kolmogorov-Smirnov distribution.
508 The two-sided Kolmogorov-Smirnov distribution has as its CDF Pr(D_n <= x),
509 for a sample of size n drawn from a distribution with CDF F(t), where
510 D_n &= sup_t |F_n(t) - F(t)|, and
511 F_n(t) is the Empirical Cumulative Distribution Function of the sample.
513 Parameters
514 ----------
515 n : integer, array_like
516 the number of samples
517 x : float, array_like
518 The K-S statistic, float between 0 and 1
519 cdf : bool, optional
520 whether to compute the CDF(default=true) or the SF.
522 Returns
523 -------
524 cdf : ndarray
525 CDF (or SF it cdf is False) at the specified locations.
527 The return value has shape the result of numpy broadcasting n and x.
528 """
529 it = np.nditer([n, x, cdf, None],
530 op_dtypes=[None, np.float64, np.bool_, np.float64])
531 for _n, _x, _cdf, z in it:
532 if np.isnan(_n):
533 z[...] = _n
534 continue
535 if int(_n) != _n:
536 raise ValueError(f'n is not integral: {_n}')
537 z[...] = _kolmogn(int(_n), _x, cdf=_cdf)
538 result = it.operands[-1]
539 return result
542def kolmognp(n, x):
543 """Computes the PDF for the two-sided Kolmogorov-Smirnov distribution.
545 Parameters
546 ----------
547 n : integer, array_like
548 the number of samples
549 x : float, array_like
550 The K-S statistic, float between 0 and 1
552 Returns
553 -------
554 pdf : ndarray
555 The PDF at the specified locations
557 The return value has shape the result of numpy broadcasting n and x.
558 """
559 it = np.nditer([n, x, None])
560 for _n, _x, z in it:
561 if np.isnan(_n):
562 z[...] = _n
563 continue
564 if int(_n) != _n:
565 raise ValueError(f'n is not integral: {_n}')
566 z[...] = _kolmogn_p(int(_n), _x)
567 result = it.operands[-1]
568 return result
571def kolmogni(n, q, cdf=True):
572 """Computes the PPF(or ISF) for the two-sided Kolmogorov-Smirnov distribution.
574 Parameters
575 ----------
576 n : integer, array_like
577 the number of samples
578 q : float, array_like
579 Probabilities, float between 0 and 1
580 cdf : bool, optional
581 whether to compute the PPF(default=true) or the ISF.
583 Returns
584 -------
585 ppf : ndarray
586 PPF (or ISF if cdf is False) at the specified locations
588 The return value has shape the result of numpy broadcasting n and x.
589 """
590 it = np.nditer([n, q, cdf, None])
591 for _n, _q, _cdf, z in it:
592 if np.isnan(_n):
593 z[...] = _n
594 continue
595 if int(_n) != _n:
596 raise ValueError(f'n is not integral: {_n}')
597 _pcdf, _psf = (_q, 1-_q) if _cdf else (1-_q, _q)
598 z[...] = _kolmogni(int(_n), _pcdf, _psf)
599 result = it.operands[-1]
600 return result