function value = tfn ( h, a ) %*****************************************************************************80 % %% TFN calculates the T function of Owen. % % Discussion: % % Owen's T function is useful for computation of the bivariate normal % distribution and the distribution of a skewed normal distribution. % % Although it was originally formulated in terms of the bivariate % normal function, the function can be defined more directly as % % T(H,A) = 1 / ( 2 * pi ) * % Integral ( 0 <= X <= A ) e^( -H^2 * (1+X^2) / 2) / (1+X^2) dX % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 26 August 2006 % % Author: % % MATLAB version by John Burkardt % % Reference: % % D B Owen, % Tables for computing the bivariate normal distribution, % Annals of Mathematical Statistics, % Volume 27, pages 1075-1090, 1956. % % J C Young and C E 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. % % Parameters: % % Input, real H, A, the arguments of the T function. % % Output, real VALUE, the value of the T function. % ngauss = 10; two_pi_inverse = 0.1591549430918953; tv1 = 1.0E-35; tv2 = 15.0; tv3 = 15.0; tv4 = 1.0E-05; weight = [ ... 0.666713443086881375935688098933E-01, ... 0.149451349150580593145776339658E+00, ... 0.219086362515982043995534934228E+00, ... 0.269266719309996355091226921569E+00, ... 0.295524224714752870173892994651E+00, ... 0.295524224714752870173892994651E+00, ... 0.269266719309996355091226921569E+00, ... 0.219086362515982043995534934228E+00, ... 0.149451349150580593145776339658E+00, ... 0.666713443086881375935688098933E-01]; xtab = [ ... -0.973906528517171720077964012084E+00, ... -0.865063366688984510732096688423E+00, ... -0.679409568299024406234327365115E+00, ... -0.433395394129247190799265943166E+00, ... -0.148874338981631210884826001130E+00, ... 0.148874338981631210884826001130E+00, ... 0.433395394129247190799265943166E+00, ... 0.679409568299024406234327365115E+00, ... 0.865063366688984510732096688423E+00, ... 0.973906528517171720077964012084E+00 ]; % % Test for H near zero. % if ( abs ( h ) < tv1 ) value = atan ( a ) * two_pi_inverse; % % Test for large values of abs(H). % elseif ( tv2 < abs ( h ) ) value = 0.0; % % Test for A near zero. % elseif ( abs ( a ) < tv1 ) value = 0.0; % % Test whether abs(A) is so large that it must be truncated. % If so, the truncated value of A is H2. % else hs = - 0.5 * h * h; h2 = a; as = a * a; % % Computation of truncation point by Newton iteration. % if ( tv3 <= log ( 1.0 + as ) - hs * as ) h1 = 0.5 * a; as = 0.25 * as; while ( 1 ) rt = as + 1.0; h2 = h1 + ( hs * as + tv3 - log ( rt ) ) ... / ( 2.0 * h1 * ( 1.0 / rt - hs ) ); as = h2 * h2; if ( abs ( h2 - h1 ) < tv4 ) break end h1 = h2; end end % % Gaussian quadrature on the interval [0,H2]. % rt = 0.0; for i = 1 : ngauss x = 0.5 * h2 * ( xtab(i) + 1.0 ); rt = rt + weight(i) * exp ( hs * ( 1.0 + x * x ) ) / ( 1.0 + x * x ); end value = rt * ( 0.5 * h2 ) * two_pi_inverse; end return end