Conflicting behaviour in Copula distribution example

I am just exploring the Copula example provided in the documentation. I noticed some strange behaviour while computing the joint pdfs. Most of the code is copied from the example without any variation - just made changes to the marginals; I input Gaussian marginals with modified means and variances as below. Both copula.prob(coordinates[0][27]) and probs[0][0][27] should ideally give the same output (which happens in most of the cases). However, sometimes (as illustrated below), copula.prob outputs inf and probs outputs nan

from __future__ import print_function
import tensorflow as tf
import tensorflow_probability as tfp
tfd = tfp.distributions

import numpy as np
import matplotlib.pyplot as plt
import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()
tfd = tfp.distributions
tfb = tfp.bijectors

class GaussianCopulaTriL(tfd.TransformedDistribution):
  """Takes a location, and lower triangular matrix for the Cholesky factor."""
  def __init__(self, loc, scale_tril):
    super(GaussianCopulaTriL, self).__init__(
        distribution=tfd.MultivariateNormalTriL(
            loc=loc,
            scale_tril=scale_tril),
        bijector=tfb.NormalCDF(),
        validate_args=False,
        name="GaussianCopulaTriLUniform")


class WarpedGaussianCopula(tfd.TransformedDistribution):
  
  def __init__(self, loc, scale_tril, marginal_bijectors, block_sizes=None):
    super(WarpedGaussianCopula, self).__init__(
        distribution=GaussianCopulaTriL(loc=loc, scale_tril=scale_tril),
        bijector=tfb.Blockwise(bijectors=marginal_bijectors,
                               block_sizes=block_sizes),
        validate_args=False,
        name="GaussianCopula")
x_axis_interval = np.linspace(-3., 3., num=30, dtype=np.float32)
y_axis_interval = np.linspace(-3., 3., num=30, dtype=np.float32)
x_grid, y_grid = np.meshgrid(x_axis_interval, y_axis_interval)

coordinates = np.concatenate(
    [x_grid[..., np.newaxis], y_grid[..., np.newaxis]], -1)


def create_gaussian_copula(correlation):
  # Use Gaussian Copula to add dependence.
  return WarpedGaussianCopula(
      loc=[0.,  0.],
      scale_tril=[[1., 0.], [correlation, tf.sqrt(1. - correlation ** 2)]],
   
      marginal_bijectors=[
          tfb.Invert(tfb.NormalCDF() (tfb.Shift(-0.87498456*0.7276069)(tfb.Scale(1/0.7276069)))  ),
          tfb.Invert(tfb.NormalCDF() (tfb.Shift(-0.0044924864*0.21490721)(tfb.Scale(1/0.2149)))) ])


# Note that the zero case will correspond to independent marginals!
correlations = [0.]
copulas = []
probs = []
for correlation in correlations:
  copula = create_gaussian_copula(correlation)
  copulas.append(copula)
  probs.append(copula.prob(coordinates))


# Plot it's density

for correlation, copula_prob in zip(correlations, probs):
  plt.figure()
  plt.contour(x_grid, y_grid, copula_prob, 100, cmap=plt.cm.jet)
  plt.title('Correlation {}'.format(correlation))

print("direct prob computation",copula.prob(coordinates[0][27])) # gives inf
print(" appended prob",probs[0][0][27]) # gives nan

No repository file or test is identified; start by running the pasted Copula example and comparing the batched copula.prob(coordinates) call with the direct copula.prob(coordinates[0][27]) call. Trace the TransformedDistribution, NormalCDF, and marginal bijector paths involved in both computations. Done means identifying why the two results differ and adding a regression test or documented reproduction for the corrected behavior.