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

96.97

/src/owens_t.jl

"""
    owens_t(x, a)

Owen's T function

Ported from Matlab implementation of John Burkardt, which is distributed under the GNU LGPL license.

Reference:
    MA Porter, DJ Winstanley, Remark AS R30: A Remark on Algorithm AS76: An Integral Useful in 
    Calculating Noncentral T and Bivariate Normal Probabilities, Applied Statistics, Volume 28, 
    Number 1, 1979, page 113.

    JC Young, Christoph Minder, Algorithm AS 76: An Algorithm Useful in Calculating Non-Central T and
    Bivariate Normal Distributions, Applied Statistics, Volume 23, Number 3, 1974, pages 455-457.
"""
function owens_t(x::Real, a::Real)
    ng = 5
    r = [0.1477621, 0.1346334, 0.1095432, 0.0747257, 0.0333357]
    tp = 0.159155
    tv1 = 1.0E-35
    tv2 = 15.0
    tv3 = 15.0
    tv4 = 1.0E-05
    u = [0.0744372, 0.2166977, 0.3397048, 0.4325317, 0.4869533]
    # Test for X near zero.
    if abs(x) < tv1
        return tp * atan(a)
    end
    # Test for large values of abs(X).
    if tv2 < abs(x)
        return 0.
    end
    # Test for `a` near zero.
    if abs(a) < tv1
        return 0.
    end
    # Test whether abs(`a`) is so large that it must be truncated.
    xs = -0.5 * x * x
    x2 = a
    fxs = a * a
    #  Computation of truncation point by Newton iteration.  
    if tv3 <= log(1.0 + fxs) - xs * fxs
        x1 = 0.5 * a
        fxs = 0.25 * fxs
        while true
            rt = fxs + 1.
            x2 = x1 + (xs * fxs + tv3 - log(rt)) / (2. * x1 * ( 1. / rt - xs))
            fxs = x2 * x2
            if abs(x2 - x1) < tv4
                break
            end
            x1 = x2
        end
    end
    # Gaussian quadrature.
    rt = 0.0
    for i in 1:ng
        r1 = 1.0 + fxs * (0.5 + u[i])^2
        r2 = 1.0 + fxs * (0.5 - u[i])^2
        rt = rt + r[i] * (exp(xs * r1) / r1 + exp(xs * r2) / r2)
    end
    return rt * x2 * tp
end

1	"""
2	owens_t(x, a)
3
4	Owen's T function
5
6	Ported from Matlab implementation of John Burkardt, which is distributed under the GNU LGPL license.
7
8	Reference:
9	MA Porter, DJ Winstanley, Remark AS R30: A Remark on Algorithm AS76: An Integral Useful in
10	Calculating Noncentral T and Bivariate Normal Probabilities, Applied Statistics, Volume 28,
11	Number 1, 1979, page 113.
12
13	JC Young, Christoph Minder, Algorithm AS 76: An Algorithm Useful in Calculating Non-Central T and
14	Bivariate Normal Distributions, Applied Statistics, Volume 23, Number 3, 1974, pages 455-457.
15	"""
16	function owens_t(x::Real, a::Real)
17	ng = 5	60,024×
18	r = [0.1477621, 0.1346334, 0.1095432, 0.0747257, 0.0333357]	300,120×
19	tp = 0.159155	60,024×
20	tv1 = 1.0E-35	60,024×
21	tv2 = 15.0	60,024×
22	tv3 = 15.0	60,024×
23	tv4 = 1.0E-05	60,024×
24	u = [0.0744372, 0.2166977, 0.3397048, 0.4325317, 0.4869533]	300,120×
25	# Test for X near zero.
26	if abs(x) < tv1	60,024×
27	return tp * atan(a)	1,224×
28	end
29	# Test for large values of abs(X).
30	if tv2 < abs(x)	58,800×
NEW 31	return 0.	!
32	end
33	# Test for `a` near zero.
34	if abs(a) < tv1	58,800×
35	return 0.	1,200×
36	end
37	# Test whether abs(`a`) is so large that it must be truncated.
38	xs = -0.5 * x * x	57,600×
39	x2 = a	57,600×
40	fxs = a * a	57,600×
41	# Computation of truncation point by Newton iteration.
42	if tv3 <= log(1.0 + fxs) - xs * fxs	57,600×
43	x1 = 0.5 * a	7,168×
44	fxs = 0.25 * fxs	7,168×
45	while true	33,024×
46	rt = fxs + 1.	33,024×
47	x2 = x1 + (xs * fxs + tv3 - log(rt)) / (2. * x1 * ( 1. / rt - xs))	33,024×
48	fxs = x2 * x2	33,024×
49	if abs(x2 - x1) < tv4	33,024×
50	break	7,168×
51	end
52	x1 = x2	25,856×
53	end
54	end
55	# Gaussian quadrature.
56	rt = 0.0	57,600×
57	for i in 1:ng	57,600×
58	r1 = 1.0 + fxs * (0.5 + u[i])^2	288,000×
59	r2 = 1.0 + fxs * (0.5 - u[i])^2	288,000×
60	rt = rt + r[i] * (exp(xs * r1) / r1 + exp(xs * r2) / r2)	518,400×
61	end
62	return rt * x2 * tp	57,600×
63	end

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