tensorflow_probability/python/internal/special_math.py
def _ndtri(p):
  """Implements ndtri core logic."""

  # Constants used in piece-wise rational approximations. Taken from the cephes
  # library:
  # https://root.cern.ch/doc/v608/SpecFuncCephesInv_8cxx_source.html
  p0 = list(reversed([-5.99633501014107895267E1,
                      9.80010754185999661536E1,
                      -5.66762857469070293439E1,
                      1.39312609387279679503E1,
                      -1.23916583867381258016E0]))
  q0 = list(reversed([1.0,
                      1.95448858338141759834E0,
                      4.67627912898881538453E0,
                      8.63602421390890590575E1,
                      -2.25462687854119370527E2,
                      2.00260212380060660359E2,
                      -8.20372256168333339912E1,
                      1.59056225126211695515E1,
                      -1.18331621121330003142E0]))
  p1 = list(reversed([4.05544892305962419923E0,
                      3.15251094599893866154E1,
                      5.71628192246421288162E1,
                      4.40805073893200834700E1,
                      1.46849561928858024014E1,
                      2.18663306850790267539E0,
                      -1.40256079171354495875E-1,
                      -3.50424626827848203418E-2,
                      -8.57456785154685413611E-4]))
  q1 = list(reversed([1.0,
                      1.57799883256466749731E1,
                      4.53907635128879210584E1,
                      4.13172038254672030440E1,
                      1.50425385692907503408E1,
                      2.50464946208309415979E0,
                      -1.42182922854787788574E-1,
                      -3.80806407691578277194E-2,
                      -9.33259480895457427372E-4]))
  p2 = list(reversed([3.23774891776946035970E0,
                      6.91522889068984211695E0,
                      3.93881025292474443415E0,
                      1.33303460815807542389E0,
                      2.01485389549179081538E-1,
                      1.23716634817820021358E-2,
                      3.01581553508235416007E-4,
                      2.65806974686737550832E-6,
                      6.23974539184983293730E-9]))
  q2 = list(reversed([1.0,
                      6.02427039364742014255E0,
                      3.67983563856160859403E0,
                      1.37702099489081330271E0,
                      2.16236993594496635890E-1,
                      1.34204006088543189037E-2,
                      3.28014464682127739104E-4,
                      2.89247864745380683936E-6,
                      6.79019408009981274425E-9]))

  def _create_polynomial(var, coeffs):
    """Compute n_th order polynomial via Horner's method."""
    coeffs = np.array(coeffs, dtype_util.as_numpy_dtype(var.dtype))
    if not coeffs.size:
      return tf.zeros_like(var)
    return coeffs[0] + _create_polynomial(var, coeffs[1:]) * var

  maybe_complement_p = tf.where(p > -np.expm1(-2.), 1. - p, p)
  # Write in an arbitrary value in place of 0 for p since 0 will cause NaNs
  # later on. The result from the computation when p == 0 is not used so any
  # number that doesn't result in NaNs is fine.
  sanitized_mcp = tf.where(
      maybe_complement_p <= 0.,
      tf.fill(tf.shape(input=p),
              dtype_util.as_numpy_dtype(p.dtype)(0.5)), maybe_complement_p)

  # Compute x for p > exp(-2): x/sqrt(2pi) = w + w**3 P0(w**2)/Q0(w**2).
  w = sanitized_mcp - 0.5
  ww = w ** 2
  x_for_big_p = w + w * ww * (_create_polynomial(ww, p0)
                              / _create_polynomial(ww, q0))
  x_for_big_p *= -np.sqrt(2. * np.pi)

  # Compute x for p <= exp(-2): x = z - log(z)/z - (1/z) P(1/z) / Q(1/z),
  # where z = sqrt(-2. * log(p)), and P/Q are chosen between two different
  # arrays based on whether p < exp(-32).
  z = tf.sqrt(-2. * tf.math.log(sanitized_mcp))
  first_term = z - tf.math.log(z) / z
  second_term_small_p = (
      _create_polynomial(1. / z, p2) /
      _create_polynomial(1. / z, q2) / z)
  second_term_otherwise = (
      _create_polynomial(1. / z, p1) /
      _create_polynomial(1. / z, q1) / z)
  x_for_small_p = first_term - second_term_small_p
  x_otherwise = first_term - second_term_otherwise

  x = tf.where(sanitized_mcp > np.exp(-2.), x_for_big_p,
               tf.where(z >= 8.0, x_for_small_p, x_otherwise))

  x = tf.where(p > 1. - np.exp(-2.), x, -x)
  infinity_scalar = tf.constant(np.inf, dtype=p.dtype)
  infinity = tf.fill(tf.shape(input=p), infinity_scalar)
  x_nan_replaced = tf.where(
      p <= 0.0, -infinity, tf.where(p >= 1.0, infinity, x))
  return x_nan_replaced
def _ndtri(p):
  """Implements ndtri core logic."""

  # Constants used in piece-wise rational approximations. Taken from the cephes
  # library:
  # https://root.cern.ch/doc/v608/SpecFuncCephesInv_8cxx_source.html
  p0 = list(reversed([-5.99633501014107895267E1,
                      9.80010754185999661536E1,
                      -5.66762857469070293439E1,
                      1.39312609387279679503E1,
                      -1.23916583867381258016E0]))
  q0 = list(reversed([1.0,
                      1.95448858338141759834E0,
                      4.67627912898881538453E0,
                      8.63602421390890590575E1,
                      -2.25462687854119370527E2,
                      2.00260212380060660359E2,
                      -8.20372256168333339912E1,
                      1.59056225126211695515E1,
                      -1.18331621121330003142E0]))
  p1 = list(reversed([4.05544892305962419923E0,
                      3.15251094599893866154E1,
                      5.71628192246421288162E1,
                      4.40805073893200834700E1,
                      1.46849561928858024014E1,
                      2.18663306850790267539E0,
                      -1.40256079171354495875E-1,
                      -3.50424626827848203418E-2,
                      -8.57456785154685413611E-4]))
  q1 = list(reversed([1.0,
                      1.57799883256466749731E1,
                      4.53907635128879210584E1,
                      4.13172038254672030440E1,
                      1.50425385692907503408E1,
                      2.50464946208309415979E0,
                      -1.42182922854787788574E-1,
                      -3.80806407691578277194E-2,
                      -9.33259480895457427372E-4]))
  p2 = list(reversed([3.23774891776946035970E0,
                      6.91522889068984211695E0,
                      3.93881025292474443415E0,
                      1.33303460815807542389E0,
                      2.01485389549179081538E-1,
                      1.23716634817820021358E-2,
                      3.01581553508235416007E-4,
                      2.65806974686737550832E-6,
                      6.23974539184983293730E-9]))
  q2 = list(reversed([1.0,
                      6.02427039364742014255E0,
                      3.67983563856160859403E0,
                      1.37702099489081330271E0,
                      2.16236993594496635890E-1,
                      1.34204006088543189037E-2,
                      3.28014464682127739104E-4,
                      2.89247864745380683936E-6,
                      6.79019408009981274425E-9]))

  def _create_polynomial(var, coeffs):
    """Compute n_th order polynomial via Horner's method."""
    coeffs = np.array(coeffs, dtype_util.as_numpy_dtype(var.dtype))
    if not coeffs.size:
      return tf.zeros_like(var)
    return coeffs[0] + _create_polynomial(var, coeffs[1:]) * var

  maybe_complement_p = tf.where(p > -np.expm1(-2.), 1. - p, p)
  # Write in an arbitrary value in place of 0 for p since 0 will cause NaNs
  # later on. The result from the computation when p == 0 is not used so any
  # number that doesn't result in NaNs is fine.
  sanitized_mcp = tf.where(
      maybe_complement_p <= 0.,
      tf.fill(tf.shape(input=p),
              dtype_util.as_numpy_dtype(p.dtype)(0.5)), maybe_complement_p)

  # Compute x for p > exp(-2): x/sqrt(2pi) = w + w**3 P0(w**2)/Q0(w**2).
  w = sanitized_mcp - 0.5
  ww = w ** 2
  x_for_big_p = w + w * ww * (_create_polynomial(ww, p0)
                              / _create_polynomial(ww, q0))
  x_for_big_p *= -np.sqrt(2. * np.pi)

  # Compute x for p <= exp(-2): x = z - log(z)/z - (1/z) P(1/z) / Q(1/z),
  # where z = sqrt(-2. * log(p)), and P/Q are chosen between two different
  # arrays based on whether p < exp(-32).
  z = tf.sqrt(-2. * tf.math.log(sanitized_mcp))
  first_term = z - tf.math.log(z) / z
  second_term_small_p = (
      _create_polynomial(1. / z, p2) /
      _create_polynomial(1. / z, q2) / z)
  second_term_otherwise = (
      _create_polynomial(1. / z, p1) /
      _create_polynomial(1. / z, q1) / z)
  x_for_small_p = first_term - second_term_small_p
  x_otherwise = first_term - second_term_otherwise

  x = tf.where(sanitized_mcp > np.exp(-2.), x_for_big_p,
               tf.where(z >= 8.0, x_for_small_p, x_otherwise))

  x = tf.where(p > 1. - np.exp(-2.), x, -x)
  infinity_scalar = tf.constant(np.inf, dtype=p.dtype)
  infinity = tf.fill(tf.shape(input=p), infinity_scalar)
  x_nan_replaced = tf.where(
      p <= 0.0, -infinity, tf.where(p >= 1.0, infinity, x))
  return x_nan_replaced
def _ndtri ( p ) : # Constants used in piece-wise rational approximations. Taken from the cephes # library: # https://root.cern.ch/doc/v608/SpecFuncCephesInv_8cxx_source.html p0 = list ( reversed ( [ - 5.99633501014107895267E1 , 9.80010754185999661536E1 , - 5.66762857469070293439E1 , 1.39312609387279679503E1 , - 1.23916583867381258016E0 ] ) ) q0 = list ( reversed ( [ 1.0 , 1.95448858338141759834E0 , 4.67627912898881538453E0 , 8.63602421390890590575E1 , - 2.25462687854119370527E2 , 2.00260212380060660359E2 , - 8.20372256168333339912E1 , 1.59056225126211695515E1 , - 1.18331621121330003142E0 ] ) ) p1 = list ( reversed ( [ 4.05544892305962419923E0 , 3.15251094599893866154E1 , 5.71628192246421288162E1 , 4.40805073893200834700E1 , 1.46849561928858024014E1 , 2.18663306850790267539E0 , - 1.40256079171354495875E-1 , - 3.50424626827848203418E-2 , - 8.57456785154685413611E-4 ] ) ) q1 = list ( reversed ( [ 1.0 , 1.57799883256466749731E1 , 4.53907635128879210584E1 , 4.13172038254672030440E1 , 1.50425385692907503408E1 , 2.50464946208309415979E0 , - 1.42182922854787788574E-1 , - 3.80806407691578277194E-2 , - 9.33259480895457427372E-4 ] ) ) p2 = list ( reversed ( [ 3.23774891776946035970E0 , 6.91522889068984211695E0 , 3.93881025292474443415E0 , 1.33303460815807542389E0 , 2.01485389549179081538E-1 , 1.23716634817820021358E-2 , 3.01581553508235416007E-4 , 2.65806974686737550832E-6 , 6.23974539184983293730E-9 ] ) ) q2 = list ( reversed ( [ 1.0 , 6.02427039364742014255E0 , 3.67983563856160859403E0 , 1.37702099489081330271E0 , 2.16236993594496635890E-1 , 1.34204006088543189037E-2 , 3.28014464682127739104E-4 , 2.89247864745380683936E-6 , 6.79019408009981274425E-9 ] ) ) def _create_polynomial ( var , coeffs ) : """Compute n_th order polynomial via Horner's method.""" coeffs = np . array ( coeffs , dtype_util . as_numpy_dtype ( var . dtype ) ) if not coeffs . size : return tf . zeros_like ( var ) return coeffs [ 0 ] + _create_polynomial ( var , coeffs [ 1 : ] ) * var maybe_complement_p = tf . where ( p > - np . expm1 ( - 2. ) , 1. - p , p ) # Write in an arbitrary value in place of 0 for p since 0 will cause NaNs # later on. The result from the computation when p == 0 is not used so any # number that doesn't result in NaNs is fine. sanitized_mcp = tf . where ( maybe_complement_p <= 0. , tf . fill ( tf . shape ( input = p ) , dtype_util . as_numpy_dtype ( p . dtype ) ( 0.5 ) ) , maybe_complement_p ) # Compute x for p > exp(-2): x/sqrt(2pi) = w + w**3 P0(w**2)/Q0(w**2). w = sanitized_mcp - 0.5 ww = w ** 2 x_for_big_p = w + w * ww * ( _create_polynomial ( ww , p0 ) / _create_polynomial ( ww , q0 ) ) x_for_big_p *= - np . sqrt ( 2. * np . pi ) # Compute x for p <= exp(-2): x = z - log(z)/z - (1/z) P(1/z) / Q(1/z), # where z = sqrt(-2. * log(p)), and P/Q are chosen between two different # arrays based on whether p < exp(-32). z = tf . sqrt ( - 2. * tf . math . log ( sanitized_mcp ) ) first_term = z - tf . math . log ( z ) / z second_term_small_p = ( _create_polynomial ( 1. / z , p2 ) / _create_polynomial ( 1. / z , q2 ) / z ) second_term_otherwise = ( _create_polynomial ( 1. / z , p1 ) / _create_polynomial ( 1. / z , q1 ) / z ) x_for_small_p = first_term - second_term_small_p x_otherwise = first_term - second_term_otherwise x = tf . where ( sanitized_mcp > np . exp ( - 2. ) , x_for_big_p , tf . where ( z >= 8.0 , x_for_small_p , x_otherwise ) ) x = tf . where ( p > 1. - np . exp ( - 2. ) , x , - x ) infinity_scalar = tf . constant ( np . inf , dtype = p . dtype ) infinity = tf . fill ( tf . shape ( input = p ) , infinity_scalar ) x_nan_replaced = tf . where ( p <= 0.0 , - infinity , tf . where ( p >= 1.0 , infinity , x ) ) return x_nan_replaced
Implements ndtri core logic.
Implements ndtri core logic .