code: def prob_mv_grid(bins, cdf, axis=-1):
      """helper function for probability of a rectangle grid in a multivariate distribution
  
      how does this generalize to more than 2 variates ?
  
      bins : tuple
          tuple of bin edges, currently it is assumed that they broadcast
          correctly
  
      """
      if not isinstance(bins, np.ndarray):
          bins = lmap(np.asarray, bins)
          n_dim = len(bins)
          bins_ = []
          # broadcast if binedges are 1d
          if all(lmap(np.ndim, bins) == np.ones(n_dim)):
              for d in range(n_dim):
                  sl = [None] * n_dim
                  sl[d] = slice(None)
                  bins_.append(bins[d][sl])
      else:  # assume it is already correctly broadcasted
          n_dim = bins.shape[0]
          bins_ = bins
  
      print(len(bins))
      cdf_values = cdf(bins_)
      probs = cdf_values.copy()
      for d in range(n_dim):
          probs = np.diff(probs, axis=d)
  
      return probs
  
imports: from statsmodels.compat.python import lmap
  import numpy
label_desc: Use image node_id 2 for calling a global function with example usage: prob_mv_grid(bins, cdf, axis) and returns: probs