JuliaMath / SpecialFunctions.jl / 814

Committed 25 Sep 2020 - 15:06 coverage: 85.701% (+0.05%) from 85.647%

Build # 814

Build Type

Pull #242

travis-ci

Committed by

web-flow

Commit Message

Merge 225521a95 into 1bb0913fc

Pull Request Pull Request #242: Adds Owen's T function

Run Details

32 of 33 new or added lines in 1 file covered. (96.97%)

62 existing lines in 10 files now uncovered.

3644 of 4252 relevant lines covered (85.7%)

2014.15 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.87

/src/beta_inc.jl

1	using Base.Math: @horner
2	using Base.MPFR: ROUNDING_MODE
3	const exparg_n = log(nextfloat(floatmin(Float64)))
4	const exparg_p = log(prevfloat(floatmax(Float64)))
5
6	#COMPUTE log(gamma(b)/gamma(a+b)) when b >= 8
7	"""
8	loggammadiv(a,b)
9
10	Computes ``log(\\Gamma(b)/\\Gamma(a+b))`` when b >= 8
11	"""
12	function loggammadiv(a::Float64, b::Float64)
13	if a > b	1,088×
14	h = b/a	!
15	c = 1.0/(1.0 + h)	!
16	x = h/(1.0 + h)	!
17	d = a + (b - 0.5)	!
18	else
19	h = a/b	1,088×
20	c = h/(1.0 + h)	1,088×
21	x = 1.0/(1.0 + h)	1,088×
22	d = b + a - 0.5	1,088×
23	end
24	x² = x*x	1,088×
25	s₃ = 1.0 + (x + x²)	1,088×
26	s₅ = 1.0 + (x + x²*s₃)	1,088×
27	s₇ = 1.0 + (x + x²*s₅)	1,088×
28	s₉ = 1.0 + (x + x²*s₇)	1,088×
29	s₁₁ = 1.0 + (x + x²*s₉)	1,088×
30
31	# SET W = stirling(b) - stirling(a+b)
32	t = inv(b)^2	1,088×
33	w = @horner(t, .833333333333333E-01, -.277777777760991E-02s₃, .793650666825390E-03s₅, -.595202931351870E-03s₇, .837308034031215E-03s₉, -.165322962780713E-02*s₁₁)	1,088×
34	w *= c/b	1,088×
35
36	#COMBINING
37	u = d*log1p(a/b)	1,088×
38	v = a*(log(b) - 1.0)	1,088×
39	return u <= v ? w - u - v : w - v - u	1,088×
40	end
41
42	"""
43	stirling_corr(a0,b0)
44
45	Compute stirling(a0) + stirling(b0) - stirling(a0 + b0)
46	for a0,b0 >= 8
47	"""
48	function stirling_corr(a0::Float64, b0::Float64)
49	a = min(a0,b0)	532×
50	b = max(a0,b0)	532×
51
52	h = a/b	532×
53	c = h/(1.0 + h)	532×
54	x = 1.0/(1.0 + h)	532×
55	x² = x*x	532×
56	#SET SN = (1-X^N)/(1-X)
57	s₃ = 1.0 + (x + x²)	532×
58	s₅ = 1.0 + (x + x²*s₃)	532×
59	s₇ = 1.0 + (x + x²*s₅)	532×
60	s₉ = 1.0 + (x + x²*s₇)	532×
61	s₁₁ = 1.0 + (x + x²*s₉)	532×
62	t = inv(b)^2	532×
63	w = @horner(t, .833333333333333E-01, -.277777777760991E-02s₃, .793650666825390E-03s₅, -.595202931351870E-03s₇, .837308034031215E-03s₉, -.165322962780713E-02*s₁₁)	532×
64	w *= c/b	532×
65	# COMPUTE stirling(a) + w
66	t = inv(a)^2	532×
67	return @horner(t, .833333333333333E-01, -.277777777760991E-02, .793650666825390E-03, -.595202931351870E-03, .837308034031215E-03, -.165322962780713E-02)/a + w	532×
68	end
69
70	"""
71	esum(mu,x)
72
73	Compute ``e^{mu+x}``
74	"""
75	function esum(mu::Float64, x::Float64)
76	if x > 0.0	1,264×
77	if mu > 0.0 \|\| mu + x < 0.0	!
78	return exp(mu)*exp(x)	!
79	else
80	return exp(mu + x)	!
81	end
82	elseif mu < 0.0 \|\| mu + x > 0.0	2,528×
83	return exp(mu)*exp(x)	200×
84	else
85	return exp(mu + x)	1,064×
86	end
87	end
88
89	"""
90	beta_integrand(a,b,x,y,mu=0.0)
91
92	Compute ``e^{mu} * x^{a}y^{b}/B(a,b)``
93	"""
94	function beta_integrand(a::Float64,b::Float64,x::Float64,y::Float64,mu::Float64=0.0)
95	a0, b0 = minmax(a,b)	1,812×
96	if a0 >= 8.0	1,264×
97	if a > b	516×
98	h = b/a	72×
99	x0 = 1.0/(1.0 + h)	72×
100	y0 = h/(1.0 + h)	72×
101	lambda = (a+b)*y - b	72×
102	else
103	h = a/b	444×
104	x0 = h/(1.0 + h)	444×
105	y0 = 1.0/(1.0 + h)	444×
106	lambda = a - (a+b)*x	444×
107	end
108	e = -lambda/a	516×
109	u = abs(e) > 0.6 ? u = e - log(x/x0) : - log1pmx(e)	876×
110	e = lambda/b	516×
111	v = abs(e) > 0.6 ? e - log(y/y0) : - log1pmx(e)	876×
112	z = esum(mu, -(au + bv))	516×
113	return (1.0/sqrt(2pi))sqrt(bx0)z*exp(-stirling_corr(a,b))	516×
114	elseif x > 0.375	748×
115	if y > 0.375	548×
116	lnx = log(x)	244×
117	lny = log(y)	244×
118	else
119	lnx = log1p(-y)	304×
120	lny = log(y)	852×
121	end
122	else
123	lnx = log(x)	200×
124	lny = log1p(-x)	200×
125	end
126	z = alnx + blny	748×
127	if a0 < 1.0	748×
128	b0 = max(a,b)	384×
129	if b0 >= 8.0	384×
130	u = loggamma1p(a0) + loggammadiv(a0,b0)	348×
131	return a0*(esum(mu, z-u))	348×
132	elseif b0 > 1.0	36×
133	u = loggamma1p(a0)	36×
134	n = trunc(Int,b0 - 1.0)	36×
135	if n >= 1	36×
UNCOV 136	c = 1.0	!
137	for i = 1:n	!
138	b0 -= 1.0	!
139	c *= (b0/(a0+b0))	!
140	end
141	u += log(c)	!
142	end
143	z -= u	36×
144	b0 -= 1.0	36×
145	apb = a0 + b0	36×
146	if apb > 1.0	36×
147	u = a0 + b0 - 1.0	36×
148	t = (1.0 + rgamma1pm1(u))/apb	36×
149	else
150	t = 1.0 + rgamma1pm1(apb)	!
151	end
152	return a0(esum(mu,z))(1.0 + rgamma1pm1(b0))/t	36×
153	end
154	else
155	z -= logbeta(a,b)	364×
156	ans = esum(mu,z)	364×
157	return ans	364×
158	end
159	if ans == 0.0	!
UNCOV 160	return 0.0	!
161	end
UNCOV 162	apb = a + b	!
UNCOV 163	if apb > 1.0	!
UNCOV 164	z = (1.0 + rgamma1pm1(apb - 1.0))/apb	!
165	else
UNCOV 166	z = 1.0 + rgamma1pm1(apb)	!
167	end
UNCOV 168	c = (1.0 + rgamma1pm1(a))*(1.0 + rgamma1pm1(b))/z	!
UNCOV 169	return ans(a0c)/(1.0 + a0/b0)	!
170	end
171
172	"""
173	beta_inc_cont_fraction(a,b,x,y,lambda,epps)
174
175	Compute ``I_{x}(a,b)`` using continued fraction expansion when `a,b > 1`.
176	It is assumed that ``\\lambda = (a+b)*y - b``
177
178	External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function)
179
180	See also: [`beta_inc`](@ref)
181
182	# Implementation
183	`BFRAC(A,B,X,Y,LAMBDA,EPS)` from Didonato and Morris (1982)
184	"""
185	function beta_inc_cont_fraction(a::Float64, b::Float64, x::Float64, y::Float64, lambda::Float64, epps::Float64)
186	ans = beta_integrand(a,b,x,y)	448×
187	if ans == 0.0	448×
188	return 0.0	!
189	end
190	c = 1.0 + lambda	448×
191	c0 = b/a	448×
192	c1 = 1.0 + 1.0/a	448×
193	yp1 = y + 1.0	448×
194
195	n = 0.0	448×
196	p = 1.0	448×
197	s = a + 1.0	448×
198	an = 0.0	448×
199	bn = 1.0	448×
200	anp1 = 1.0	448×
201	bnp1 = c/c1	448×
202	r = c1/c	448×
203	#CONT FRACTION
204
205	while true	4,356×
206	n += 1.0	4,356×
207	t = n/a	4,356×
208	w = n(b - n)x	4,356×
209	e = a/s	4,356×
210	alpha = (p(p+c0)ee)(w*x)	4,356×
211	e = (1.0 + t)/(c1 + 2*t)	4,356×
212	beta = n + w/s +e(c + nyp1)	4,356×
213	p = 1.0 + t	4,356×
214	s += 2.0	4,356×
215
216	#update an, bn, anp1, bnp1
217	t = alphaan + betaanp1	4,356×
218	an = anp1	4,356×
219	anp1 = t	4,356×
220	t = alphabn + betabnp1	4,356×
221	bn = bnp1	4,356×
222	bnp1 = t	4,356×
223
224	r0 = r	4,356×
225	r = anp1/bnp1	4,356×
226	if abs(r - r0) <= epps*r	4,356×
227	break	448×
228	end
229	#rescale
230	an /= bnp1	3,908×
231	bn /= bnp1	3,908×
232	anp1 = r	3,908×
233	bnp1 = 1.0	3,908×
234	end
235	return ans*r	448×
236	end
237
238	"""
239	beta_inc_asymptotic_symmetric(a,b,lambda,epps)
240
241	Compute ``I_{x}(a,b)`` using asymptotic expansion for `a,b >= 15`.
242	It is assumed that ``\\lambda = (a+b)*y - b``.
243
244	External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function)
245
246	See also: [`beta_inc`](@ref)
247
248	# Implemention
249	`BASYM(A,B,LAMBDA,EPS)` from Didonato and Morris (1982)
250	"""
251	function beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64, epps::Float64)
252	a0 =zeros(22)	264×
253	b0 = zeros(22)	264×
254	c = zeros(22)	264×
255	d = zeros(22)	264×
256	e0 = 2/sqrt(pi)	12×
257	e1 = 2^(-1.5)	12×
258	sm = 0.0	12×
259	ans = 0.0	12×
260	if a > b	12×
261	h = b/a	!
262	r0 = 1.0/(1.0 + h)	!
263	r1 = (b-a)/a	!
264	w0 = 1.0/sqrt(b*(1.0+h))	!
265	else
266	h = a/b	12×
267	r0 = 1.0/(1.0 + h)	12×
268	r1 = (b-a)/b	12×
269	w0 = 1.0/sqrt(a*(1.0+h))	12×
270	end
271	f = -alog1pmx(-(lambda/a)) - blog1pmx((lambda/b))	12×
272	t = exp(-f)	12×
273	if t == 0.0	12×
274	return ans	!
275	end
276	z0 = sqrt(f)	12×
277	z = 0.5*(z0/e1)	12×
278	z² = 2.0*f	12×
279
280	a0[1] = (2.0/3.0)*r1	12×
281	c[1] = -0.5*a0[1]	12×
282	d[1] = - c[1]	12×
283	j0 = (0.5/e0)*erfcx(z0)	12×
284	j1 = e1	12×
285	sm = j0 + d[1]w0j1	12×
286
287	s = 1.0	12×
288	h² = h*h	12×
289	hn = 1.0	12×
290	w = w0	12×
291	znm1 = z	12×
292	zn = z²	12×
293
294	for n = 2: 2: 20	12×
295	hn *= h²	72×
296	a0[n] = 2.0r0(1.0 + h*hn)/(n + 2.0)	72×
297	s += hn	72×
298	a0[n+1] = 2.0r1s/(n+3.0)	72×
299
300	for i = n: n+1	144×
301	r = -0.5*(i + 1.0)	144×
302	b0[1] = r*a0[1]	144×
303	for m = 2:i	288×
304	bsum = 0.0	936×
305	for j =1: m-1	1,872×
306	bsum += (jr - (m-j))a0[j]*b0[m-j]	7,800×
307	end
308	b0[m] = r*a0[m] + bsum/m	1,728×
309	end
310	c[i] = b0[i]/(i+1.0)	144×
311	dsum = 0.0	144×
312	for j = 1: i-1	288×
313	imj = i - j	936×
314	dsum += d[imj]*c[j]	1,728×
315	end
316	d[i] = -(dsum + c[i])	216×
317	end
318
319	j0 = e1znm1 + (n - 1)j0	72×
320	j1 = e1zn + nj1	72×
321	znm1 *= z²	72×
322	zn *= z²	72×
323	w *= w0	72×
324	t0 = d[n]wj0	72×
325	w *= w0	72×
326	t1 = d[n+1]wj1	72×
327	sm += (t0 + t1)	72×
328	if (abs(t0) + abs(t1)) <= epps*sm	72×
329	break	72×
330	end
331	end
332
333	u = exp(-stirling_corr(a,b))	12×
334	return e0tu*sm	12×
335	end
336
337	"""
338	beta_inc_asymptotic_asymmetric(a,b,x,y,w,epps)
339
340	Evaluation of ``I_{x}(a,b)`` when `b < min(epps,epps*a)` and `x <= 0.5` using asymptotic expansion.
341	It is assumed `a >= 15` and `b <= 1`, and epps is tolerance used.
342
343	External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function)
344
345	See also: [`beta_inc`](@ref)
346	"""
347	function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::Float64, w::Float64, epps::Float64)
348	c = zeros(31)	10,292×
349	d = zeros(31)	10,292×
350	bm1 = b - 1.0	332×
351	nu = a + 0.5*bm1	332×
352	if y > 0.375	332×
353	lnx = log(x)	!
354	else
355	lnx = log1p(-y)	332×
356	end
357	z = -nu*lnx	332×
358	if b*z == 0.0	332×
359	return error("expansion can't be computed")	!
360	end
361
362	# COMPUTATION OF THE EXPANSION
363	#SET R = EXP(-Z)Z*B/GAMMA(B)
364	r = b(1.0 + rgamma1pm1(b))exp(b*log(z))	332×
365	r = exp(alnx)exp(0.5bm1*lnx)	332×
366	u = loggammadiv(b,a) + b*log(nu)	332×
367	u = r*exp(-u)	332×
368	if u == 0.0	332×
369	return error("expansion can't be computed")	!
370	end
371	(p, q) = gamma_inc(b,z,0)	332×
372	v = inv(nu)^2/4	332×
373	t2 = lnx^2/4	332×
374	l = w/u	332×
375	j = q/r	332×
376	sm = j	332×
377	t = 1.0	332×
378	cn = 1.0	332×
379	n2 = 0.0	332×
380	for n = 1:30	332×
381	bp2n = b + n2	1,184×
382	j = (bp2n(bp2n + 1.0)j + (z + bp2n + 1.0)t)v	1,184×
383	n2 += 2.0	1,184×
384	t *= t2	1,184×
385	cn /= n2*(n2 + 1.0)	1,184×
386	c[n] = cn	1,184×
387	s = 0.0	1,184×
388	if n != 1	1,184×
389	nm1 = n -1	852×
390	coef = b - n	852×
391	for i = 1:nm1	1,704×
392	s += coefc[i]d[n-i]	2,112×
393	coef += b	3,372×
394	end
395	end
396	d[n] = bm1*cn + s/n	1,184×
397	dj = d[n] * j	1,184×
398	sm += dj	1,184×
399	if sm <= 0.0	1,184×
400	return error("expansion can't be computed")	!
401	end
402	if abs(dj) <= epps*(sm+l)	1,184×
403	break	1,184×
404	end
405	end
406	return w + u*sm	332×
407	end
408
409
410	#For b < min(eps, eps*a) and x <= 0.5
411	"""
412	beta_inc_power_series2(a,b,x,epps)
413
414	Variant of `BPSER(A,B,X,EPS)`.
415
416	External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function)
417
418	See also: [`beta_inc`](@ref)
419
420	# Implementation
421	`FPSER(A,B,X,EPS)` from Didonato and Morris (1982)
422	"""
423	function beta_inc_power_series2(a::Float64, b::Float64, x::Float64, epps::Float64)
424	ans = 1.0	!
425	if a > 1.0e-3 * epps	!
UNCOV 426	ans = 0.0	!
427	t = a*log(x)	!
428	if t < exparg_n	!
429	return ans	!
430	end
431	ans = exp(t)	!
432	end
433	ans *= b/a	!
434	tol = epps/a	!
435	an = a + 1.0	!
UNCOV 436	t = x	!
437	s = t/an	!
438	an += 1.0	!
439	t *= x	!
440	c = t/an	!
441	s += c	!
442	while abs(c) > tol	!
443	an += 1.0	!
444	t *= x	!
445	c = t/an	!
446	s += c	!
447	end
448	ans = (1.0 + as)	!
449	return ans	!
450	end
451
452	#A <= MIN(EPS,EPSB), BX <= 1, AND X <= 0.5., A is small
453	"""
454	beta_inc_power_series1(a,b,x,epps)
455
456	Another variant of `BPSER(A,B,X,EPS)`.
457
458	External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function)
459
460	See also: [`beta_inc`](@ref)
461
462	# Implementation
463	`APSER(A,B,X,EPS)` from Didonato and Morris (1982)
464	"""
465	function beta_inc_power_series1(a::Float64, b::Float64, x::Float64, epps::Float64)
466	g = Base.MathConstants.eulergamma	!
467	bx = b*x	!
468	t = x - bx	!
469	if b*epps > 2e-2	!
470	c = log(bx) + g + t	!
471	else
472	c = log(x) + digamma(b) + g + t	!
473	end
474	tol = 5.0eppsabs(c)	!
UNCOV 475	j = 1.0	!
UNCOV 476	s = 0.0	!
UNCOV 477	j += 1.0	!
478	t *= (x - bx/j)	!
479	aj = t/j	!
480	s += aj	!
481	while abs(aj) > tol	!
482	j += 1.0	!
483	t *= (x - bx/j)	!
484	aj = t/j	!
485	s += aj	!
486	end
487	return -a*(c + s)	!
488	end
489
490	#B .LE. 1 OR B*X .LE. 0.7
491	"""
492	beta_inc_power_series(a,b,x,epps)
493
494	Computes ``I_x(a,b)`` using power series :
495	```math
496	I_{x}(a,b) = G(a,b)x^{a}/a (1 + a\\sum_{j=1}^{\\infty}((1-b)(2-b)...(j-b)/j!(a+j)) x^{j})
497	```
498	External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function)
499
500	See also: [`beta_inc`](@ref)
501
502	# Implementation
503	`BPSER(A,B,X,EPS)` from Didonato and Morris (1982)
504	"""
505	function beta_inc_power_series(a::Float64, b::Float64, x::Float64, epps::Float64)
506	ans = 0.0	2,356×
507	if x == 0.0	2,356×
508	return 0.0	!
509	end
510	a0 = min(a,b)	2,356×
511	b0 = max(a,b)	2,356×
512	if a0 >= 1.0	2,356×
513	z = a*log(x) - logbeta(a,b)	604×
514	ans = exp(z)/a	604×
515	else
516
517	if b0 >= 8.0	1,752×
518	u = loggamma1p(a0) + loggammadiv(a0,b0)	396×
519	z = a*log(x) - u	396×
520	ans = (a0/a)*exp(z)	396×
521	if ans == 0.0 \|\| a <= 0.1*epps	792×
522	return ans	!
523	end
524	elseif b0 > 1.0	1,356×
525	u = loggamma1p(a0)	36×
526	m = b0 - 1.0	36×
527	if m >= 1.0	36×
UNCOV 528	c = 1.0	!
529	for i = 1:m	!
530	b0 -= 1.0	!
531	c *= (b0/(a0+b0))	!
532	end
533	u += log(c)	!
534	end
535	z = a*log(x) - u	36×
536	b0 -= 1.0	36×
537	apb = a0 + b0	36×
538	if apb > 1.0	36×
539	u = a0 + b0 - 1.0	36×
540	t = (1.0 + rgamma1pm1(u))/apb	36×
541	else
542	t = 1.0 + rgamma1pm1(apb)	!
543	end
544	ans = exp(z)(a0/a)(1.0 + rgamma1pm1(b0))/t	36×
545	if ans == 0.0 \|\| a <= 0.1*epps	72×
546	return ans	!
547	end
548	else
549	#PROCEDURE FOR A0 < 1 && B0 < 1
550	ans = x^a	1,320×
551	if ans == 0.0	1,320×
552	return ans	!
553	end
554	apb = a + b	1,320×
555	if apb > 1.0	1,320×
556	u = a + b - 1.0	892×
557	z = (1.0 + rgamma1pm1(u))/apb	892×
558	else
559	z = 1.0 + rgamma1pm1(apb)	428×
560	end
561	c = (1.0 + rgamma1pm1(a))*(1.0 + rgamma1pm1(b))/z	1,320×
562	ans = c(b/apb)	1,320×
563	#label l70 start
564	if ans == 0.0 \|\| a <= 0.1*epps	2,640×
565	return ans	!
566	end
567	end
568	end
569	if ans == 0.0 \|\| a <= 0.1*epps	4,712×
570	return ans	!
571	end
572	# COMPUTE THE SERIES
573
574	sm = 0.0	2,356×
575	n = 0.0	2,356×
576	c = 1.0	2,356×
577	tol = epps/a	2,356×
578	n += 1.0	2,356×
579	c = x(1.0 - b/n)	2,356×
580	w = c/(a + n)	2,356×
581	sm += w	2,356×
582	while abs(w) > tol	48,000×
583	n += 1.0	45,644×
584	c = x(1.0 - b/n)	45,644×
585	w = c/(a+n)	45,644×
586	sm += w	45,644×
587	end
588	return ans(1.0 + asm)	2,356×
589	end
590
591	"""
592	beta_inc_diff(a,b,x,y,n,epps)
593
594	Compute ``I_{x}(a,b) - I_{x}(a+n,b)`` where `n` is positive integer and epps is tolerance.
595	A generalised version of [DLMF](https://dlmf.nist.gov/8.17.20).
596
597	External links: [DLMF](https://dlmf.nist.gov/8.17.20), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function)
598
599	See also: [`beta_inc`](@ref)
600	"""
601	function beta_inc_diff(a::Float64, b::Float64, x::Float64, y::Float64, n::Integer, epps::Float64)
602	apb = a + b	716×
603	ap1 = a + 1.0	716×
604	mu = 0.0	716×
605	d = 1.0	716×
606	if n != 1 && a >= 1.0 && apb >= 1.1*ap1	716×
607	mu = abs(exparg_n)	200×
608	k = exparg_p	200×
609	if k < mu	200×
UNCOV 610	mu = k	!
611	end
612	t = mu	200×
613	d = exp(-t)	200×
614	end
615
616	ans = beta_integrand(a,b,x,y,mu)/a	716×
617	if n == 1 \|\| ans == 0.0	1,148×
618	return ans	284×
619	end
620	nm1 = n -1	432×
621	w = d	432×
622
623	k = 0	432×
624	if b <= 1.0	432×
625	kp1 = k + 1	196×
626	for i = kp1:nm1	392×
627	l = i - 1	3,724×
628	d = ((apb + l)/(ap1 + l))x	3,724×
629	w += d	3,724×
630	if d <= epps*w	3,724×
631	break	3,724×
632	end
633	end
634	return ans*w	196×
635	elseif y > 1.0e-4	236×
636	r = trunc(Int,(b - 1.0)*x/y - a)	236×
637	if r < 1.0	236×
UNCOV 638	kp1 = k + 1	!
639	for i = kp1:nm1	!
640	l = i - 1	!
641	d = ((apb + l)/(ap1 + l))x	!
642	w += d	!
643	if d <= epps*w	!
644	break	!
645	end
646	end
647	return ans*w	!
648	end
649	k = t = nm1	236×
650	if r < t	236×
651	k = r	44×
652	end
653	# ADD INC TERMS OF SERIES
654	for i = 1:k	472×
655	l = i -1	2,508×
656	d = ((apb + l)/(ap1 + l))x	2,508×
657	w += d	4,780×
658	end
659	if k == nm1	236×
660	return ans*w	192×
661	end
662	else
UNCOV 663	k = nm1	!
664	for i = 1:k	!
665	l = i -1	!
666	d = ((apb + l)/(ap1 + l))x	!
667	w += d	!
668	end
669	if k == nm1	!
670	return ans*w	!
671	end
672	end
673
674	kp1 = k + 1	44×
675	for i = kp1:nm1	88×
676	l = i - 1	172×
677	d = ((apb + l)/(ap1 + l))x	172×
678	w += d	172×
679	if d <= epps*w	172×
680	break	172×
681	end
682	end
683	return ans*w	44×
684	end
685
686
687	#SIGNIFICANT DIGIT COMPUTATION OF INCOMPLETE BETA FUNCTION RATIOS
688	#by ARMIDO R. DIDONATO AND ALFRED H. MORRIS, JR.
689	#ACM Transactions on Mathematical Software. Vol18, No3, September1992, Pages360-373
690	#DLMF : https://dlmf.nist.gov/8.17#E1
691	#Wikipedia : https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function
692
693	"""
694	beta_inc(a,b,x)
695
696	Returns a tuple ``(I_{x}(a,b),1.0-I_{x}(a,b))`` where the Regularized Incomplete Beta Function is given by:
697	```math
698	I_{x}(a,b) = \\frac{1}{B(a,b)} \\int_{0}^{x} t^{a-1}(1-t)^{b-1} dt,
699	```
700	where ``B(a,b) = \\Gamma(a)\\Gamma(b)/\\Gamma(a+b)``.
701
702	External links: [DLMF](https://dlmf.nist.gov/8.17.1), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function)
703
704	See also: [`beta_inc_inv(a,b,p,q)`](@ref SpecialFunctions.beta_inc_inv)
705	"""
706	function beta_inc(a::Float64, b::Float64, x::Float64, y::Float64)
707	p = 0.0	3,148×
708	q = 0.0	3,148×
709	# lambda = a - (a+b)*x
710	if a < 0.0 \|\| b < 0.0	6,292×
711	return error("a or b is negative")	!
712	elseif a == 0.0 && b == 0.0	3,148×
713	return error("a and b are 0.0")	!
714	elseif x < 0.0 \|\| x > 1.0	6,292×
715	return error("x < 0 or x > 1")	!
716	elseif y < 0.0 \|\| y > 1.0	6,292×
717	return error("y < 0 or y > 1")	!
718	else
719	z = x + y - 1.0	3,148×
720	if abs(z) > 3.0*eps()	3,148×
721	return error("x + y != 1.0") # ERROR HANDLING	!
722	end
723	end
724
725	if x == 0.0	3,148×
726	return (0.0,1.0)	!
727	elseif y == 0.0	3,148×
728	return (1.0,0.0)	!
729	elseif a == 0.0	3,148×
730	return (1.0,0.0)	!
731	elseif b == 0.0	3,148×
732	return (0.0,1.0)	!
733	end
734	#EVALUATION OF ALGOS FOR PROPER SUB-DOMAINS ABOVE
735	epps = max(eps(), 1.0e-15)	3,148×
736	if max(a,b) < 1.0E-3 * epps	3,148×
737	return (b/(a+b), a/(a+b))	!
738	end
739	ind = false	3,148×
740	a0 = a	3,148×
741	b0 = b	3,148×
742	x0 = x	3,148×
743	y0 = y	3,148×
744
745	if min(a0,b0) > 1.0	3,148×
746	#PROCEDURE FOR A0>1 AND B0>1
747	lambda = a > b ? (a+b)y - b : a - (a+b)x	2,432×
748	if lambda < 0.0	1,356×
749	ind = true	604×
750	a0 = b	604×
751	b0 = a	604×
752	x0 = y	604×
753	y0 = x	604×
754	lambda = abs(lambda)	604×
755	end
756	if b0 < 40.0 && b0*x0 <= 0.7	1,356×
757	p = beta_inc_power_series(a0,b0,x0,epps)	376×
758	q = 1.0 - p	376×
759	elseif b0 < 40.0	980×
760	n = trunc(Int, b0)	520×
761	b0 -= float(n)	520×
762	if b0 == 0.0	520×
763	n-=1	276×
764	b0=1.0	276×
765	end
766	p = beta_inc_diff(b0,a0,y0,x0,n,epps)	520×
767	if x0 <= 0.7	520×
768	p += beta_inc_power_series(a0,b0,x0,epps)	372×
769	q = 1.0 - p	372×
770	else
771	if a0 <= 15.0	148×
772	n = 20	88×
773	p += beta_inc_diff(a0, b0, x0, y0, n, epps)	88×
774	a0 += n	88×
775	end
776	p = beta_inc_asymptotic_asymmetric(a0,b0,x0,y0,p,15.0*eps())	148×
777	q = 1.0 - p	664×
778	end
779	elseif a0 > b0	460×
780	if b0 <= 100.0 \|\| lambda > 0.03*b0	!
781	p = beta_inc_cont_fraction(a0,b0,x0,y0,lambda,15.0*eps())	!
782	q = 1.0 - p	!
783	else
784	p = beta_inc_asymptotic_symmetric(a0,b0,lambda,100.0*eps())	!
785	q = 1.0 - p	!
786	end
787	elseif a0 <= 100.0 \|\| lambda > 0.03*a0	856×
788	p = beta_inc_cont_fraction(a0,b0,x0,y0,lambda,15.0*eps())	448×
789	q = 1.0 - p	448×
790	else
791	p = beta_inc_asymptotic_symmetric(a0,b0,lambda,100.0*eps())	12×
792	q = 1.0 - p	12×
793	end
794	return ind ? (q, p) : (p, q)	1,356×
795	end
796	#PROCEDURE FOR A0<=1 OR B0<=1
797	if x > 0.5	1,792×
798	ind = true	832×
799	a0 = b	832×
800	b0 = a	832×
801	y0 = x	832×
802	x0 = y	832×
803	end
804
805	if b0 < min(epps, epps*a0)	1,792×
806	p = beta_inc_power_series2(a0,b0,x0,epps)	!
807	q = 1.0 - p	!
808	elseif a0 < min(epps, eppsb0) && b0x0 <= 1.0	1,792×
809	q = beta_inc_power_series1(a0,b0,x0,epps)	!
810	p = 1.0 - q	!
811	elseif max(a0,b0) > 1.0	1,792×
812	if b0 <= 1.0	464×
813	p = beta_inc_power_series(a0,b0,x0,epps)	200×
814	q = 1.0 - p	200×
815	elseif x0 >= 0.3	264×
816	q = beta_inc_power_series(b0,a0,y0,epps)	80×
817	p = 1.0 - q	80×
818	elseif x0 >= 0.1	184×
819	if b0 > 15.0	132×
820	q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())	76×
821	p = 1.0 - q	76×
822	else
823	n = 20	56×
824	q = beta_inc_diff(b0,a0,y0,x0,n,epps)	56×
825	b0 += n	56×
826	q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())	56×
827	p = 1.0 - q	188×
828	end
829	elseif (x0*b0)^(a0) <= 0.7	52×
830	p = beta_inc_power_series(a0,b0,x0,epps)	!
831	q = 1.0 - p	!
832	else
833	n = 20	52×
834	q = beta_inc_diff(b0,a0,y0,x0,n,epps)	52×
835	b0 += n	52×
836	q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())	52×
837	p = 1.0 - q	516×
838	end
839	elseif a0 >= min(0.2, b0)	1,328×
840	p = beta_inc_power_series(a0,b0,x0,epps)	1,328×
841	q = 1.0 - p	1,328×
842	elseif x0^a0 <= 0.9	!
843	p = beta_inc_power_series(a0,b0,x0,epps)	!
844	q = 1.0 - p	!
845	elseif x0 >= 0.3	!
846	q = beta_inc_power_series(b0,a0,y0,epps)	!
847	p = 1.0 - q	!
848	else
UNCOV 849	n = 20	!
850	q = beta_inc_diff(b0,a0,y0,x0,n,epps)	!
851	b0 += n	!
852	q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())	!
853	p = 1.0 - q	!
854	end
855
856	#TERMINATION
857	return ind ? (q, p) : (p, q)	1,792×
858	end
859
860	beta_inc(a::Float64, b::Float64, x::Float64) = beta_inc(a, b, x, 1.0 - x)	768×
UNCOV 861	function beta_inc(a::T, b::T, x::T, y::T) where {T<:Union{Float16, Float32}}	!
862	T.(beta_inc(Float64(a), Float64(b), Float64(x), Float64(y)))	!
863	end
864	beta_inc(a::Real, b::Real, x::Real, y::Real) = beta_inc(promote(float(a), float(b), float(x), float(y))...)	!
865	beta_inc(a::T, b::T, x::T, y::T) where {T<:AbstractFloat} = throw(MethodError(beta_inc,(a, b, x, y,"")))	!
866
867	#GW Cran, KJ Martin, GE Thomas, Remark AS R19 and Algorithm AS 109: A Remark on Algorithms AS 63: The Incomplete Beta Integral and AS 64: Inverse of the Incomplete Beta Integeral,
868	#Applied Statistics,
869	#Volume 26, Number 1, 1977, pages 111-114.
870
871	"""
872	beta_inc_inv(a,b,p,q,lb=logbeta(a,b)[1])
873
874	Computes inverse of incomplete beta function. Given `a`,`b` and ``I_{x}(a,b) = p`` find `x` and return tuple `(x,y)`.
875
876	See also: [`beta_inc(a,b,x)`](@ref SpecialFunctions.beta_inc)
877	"""
878	function beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64; lb = logbeta(a,b)[1])
879	fpu = 1e-30	360×
880	x = p	180×
881	if p == 0.0	180×
882	return (0.0, 1.0)	4×
883	elseif p == 1.0	176×
884	return (1.0, 0.0)	8×
885	end
886
887	#change tail if necessary
888
889	if p > 0.5	168×
890	pp = q	44×
891	aa = b	44×
892	bb = a	44×
893	indx = true	44×
894	else
895	pp = p	124×
896	aa = a	124×
897	bb = b	124×
898	indx = false	124×
899	end
900
901	#Initial approx
902
903	r = sqrt(-log(pp^2))	168×
904	pp_approx = r - @horner(r, 2.30753e+00, 0.27061e+00) / @horner(r, 1.0, .99229e+00, .04481e+00)	168×
905
906	if a > 1.0 && b > 1.0	168×
907	r = (pp_approx^2 - 3.0)/6.0	108×
908	s = 1.0/(2*aa - 1.0)	108×
909	t = 1.0/(2*bb - 1.0)	108×
910	h = 2.0/(s+t)	108×
911	w = pp_approxsqrt(h+r)/h - (t-s)(r + 5.0/6.0 - 2.0/(3.0*h))	108×
912	x = aa/ (aa+bb*exp(w^2))	108×
913	else
914	r = 2.0*bb	60×
915	t = 1.0/(9.0*bb)	60×
916	t = r(1.0-t+pp_approxsqrt(t))^3	60×
917	if t <= 0.0	60×
918	x = -expm1((log((1.0-pp)*bb)+lb)/bb)	!
919	else
920	t = (4.0*aa+r-2.0)/t	60×
921	if t <= 1.0	60×
922	x = exp((log(pp*aa)+lb)/aa)	20×
923	else
924	x = 1.0 - 2.0/(t+1.0)	40×
925	end
926	end
927	end
928
929	#solve x using modified newton-raphson iteration
930
931	r = 1.0 - aa	168×
932	t = 1.0 - bb	168×
933	pp_approx_prev = 0.0	168×
934	sq = 1.0	168×
935	prev = 1.0	168×
936
937	if x < 0.0001	168×
UNCOV 938	x = 0.0001	!
939	end
940	if x > .9999	168×
UNCOV 941	x = .9999	!
942	end
943
944	iex = max(-5.0/aa^2 - 1.0/pp^0.2 - 13.0, -30.0)	168×
945	acu = 10.0^iex	168×
946
947	#iterate
948	while true	668×
949	pp_approx = beta_inc(aa,bb,x)[1]	668×
950	xin = x	668×
951	pp_approx = (pp_approx-pp)exp(lb+rlog(xin)+t*log1p(-xin))	668×
952	if pp_approx * pp_approx_prev <= 0.0	668×
953	prev = max(sq, fpu)	232×
954	end
955	g = 1.0	668×
956
957	tx = x - g*pp_approx	668×
958	while true	668×
959	adj = g*pp_approx	668×
960	sq = adj^2	668×
961	tx = x - adj	668×
962	if (prev > sq && tx >= 0.0 && tx <= 1.0)	668×
963	break	668×
964	end
965	g /= 3.0	!
966	end
967
968	#check if current estimate is acceptable
969
970	if prev <= acu \|\| pp_approx^2 <= acu	1,336×
971	x = tx	168×
972	return indx ? (1.0 - x, x) : (x, 1.0-x)	168×
973	end
974
975	if tx == x	500×
976	return indx ? (1.0 - x, x) : (x, 1.0-x)	!
977	end
978
979	x = tx	500×
980	pp_approx_prev = pp_approx	500×
981	end
982	end
983
984	beta_inc_inv(a::Float64, b::Float64, p::Float64) = beta_inc_inv(a, b, p, 1.0-p)	180×

JuliaMath / SpecialFunctions.jl / 814

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous