R/tOwenT.R

In Phi2rho: Owen's T Function and Bivariate Normal Integral

tOwenT <- function( h, r, ph, prh, transf = TRUE, pi ) {
  # This function calculates the Owen''s function T( h, r ). It can handle the
  # special cases r = Inf and r = -Inf. In such cases, indefinite expressions
  # occur but are caught and modified so that the loop stops with s = 0, and
  # then the final result is formed. The calling function ensures that h and
  # r are of the same length and provides pnorm( h ) and pnorm( r * h ) in ph
  # and prh parameters respectively, where prh must be set to 0 if h is zero
  # and r is infinite.
  maxIter <- getOption( "Phi2rho.maxIter", default = 5000 )
  u <- ph 
  # i indicates where the transformations must be done
  i <- ( transf & abs( r ) > 1 ) | is.infinite( r )
  # attention: the transformation of u must be before the transformation of r,
  # because otherwise for r = Inf and r = -Inf we get the sign 0 instead of 1
  # and -1 repectively, because sign( 0 ) = 0 in R
  u[ i ] <- ( u / 2 + prh * ( 1 / 2 - u ) - ( 1 - sign( r ) ) / 4 )[ i ]
  h[ i ] <- r[ i ] * h[ i ]
  r[ i ] <- 1 / r[ i ]
  # Parameter checking in the calling function ensures that h and r have the
  # same precBits if they are mpfr numbers. If Rmpfr is used, all variables
  # below are automatically given class mpfr due to calculation with h or r,
  # or with other variables that have already inherited the mpfr class.
  He0 <- rep( 1, length( h ) ) # = He_0(h)
  He1 <- h                     # = He_1(h)
  p <- r^2 / ( 1 + r^2 )
  q <- sqrt( p ) * dnorm( h ) / sqrt( 2 * pi )
  s <- q
  z <- s * 0
  n <- rep( -1, length( s ) ) # number of iterations for each component
  k <- 0
  repeat {
    k <- k + 1
    He0 <- h * He1 - ( 2 * k - 1 ) * He0 # = He_{2k} / ( 2^k * k! )
    He1 <- h * He0 - 2 * k * He1 # = He_{2k+1} / ( 2^k * k! )
    He0 <- He0 / ( 2 * k )
    He1 <- He1 / ( 2 * k )
    q <- -q * p
    v <- s + He0 * q / ( 2 * k + 1 )
    v[ is.nan( v ) ] <- 0 # needed because of special cases
    # if h is a zero of the 2k-th Hermite polynomial, the new approximation is
    # equal to the last one, but the recursion must not be terminated, so the
    # recursion is terminated when the new approximation is equal to the last
    # one and penultimate one
    n[ n == -1 & v == s & s == z ] <- k
    if ( all( n > 0 ) || k >= maxIter ) break
    z <- s
    s <- v
  }
  s[ r < 0 ] <- -s[ r < 0 ] # equating the sign with the sign of r
  s[ i ] <- u[ i ] - s[ i ] # adjustment due to parameter transformation
  attr( s, "nIter" ) <- n
  return( s )
}