code: def multivariate_normal(x, mu, L):
      
      mu_ndims = mu.shape.ndims
      x_ndims = x.shape.ndims
      if x_ndims is not None and x_ndims != 2:
          raise ValueError('Shape of x must be 2D.')
      if mu_ndims is not None and mu_ndims != 2:
          raise ValueError('Shape of mu must be 2D.')
  
      d = x - mu
      alpha = tf.matrix_triangular_solve(L, d, lower=True)
      num_dims = tf.cast(tf.shape(d)[0], L.dtype)
      p = - 0.5 * tf.reduce_sum(tf.square(alpha), 0)
      p -= 0.5 * num_dims * np.log(2 * np.pi)
      p -= tf.reduce_sum(tf.log(tf.matrix_diag_part(L)))
      return p
docs: Computes the log-density of a multivariate normal.
      :param x  : Dx1 or DxN sample(s) for which we want the density
      :param mu : Dx1 or DxN mean(s) of the normal distribution
      :param L  : DxD Cholesky decomposition of the covariance matrix
      :return p : (1,) or (N,) vector of log densities for each of the N x's and/or mu's
  
      x and mu are either vectors or matrices. If both are vectors (N,1):
      p[0] = log pdf(x) where x ~ N(mu, LL^T)
      If at least one is a matrix, we assume independence over the *columns*:
      the number of rows must match the size of L. Broadcasting behaviour:
      p[n] = log pdf of:
      x[n] ~ N(mu, LL^T) or x ~ N(mu[n], LL^T) or x[n] ~ N(mu[n], LL^T)