code: def probability_density(self, X):
          
          self.check_fit()
  
          U, V = self.split_matrix(X)
  
          if self.theta == 1:
              return np.multiply(U, V)
  
          else:
              a = np.power(np.multiply(U, V), -1)
              tmp = np.power(-np.log(U), self.theta) + np.power(-np.log(V), self.theta)
              b = np.power(tmp, -2 + 2.0 / self.theta)
              c = np.power(np.multiply(np.log(U), np.log(V)), self.theta - 1)
              d = 1 + (self.theta - 1) * np.power(tmp, -1.0 / self.theta)
              return self.cumulative_distribution(X) * a * b * c * d
docs: Compute density function for given copula family.