code: def probability_density(self, X):
          
          self.check_fit()
  
          U, V = self.split_matrix(X)
  
          if self.theta == 0:
              return np.multiply(U, V)
  
          else:
              num = np.multiply(np.multiply(-self.theta, self._g(1)), 1 + self._g(np.add(U, V)))
              aux = np.multiply(self._g(U), self._g(V)) + self._g(1)
              den = np.power(aux, 2)
              return num / den
docs: Compute density function for given copula family.
  
          Args:
              X: `np.ndarray`
  
          Returns:
              np.array: probability density