code: def _statsmodels_bivariate_kde(self, x, y, bws, xsize, ysize, xyranges):
          
          import statsmodels.nonparametric.api as smnp
          for i, (coord, bw) in enumerate(zip([x, y], bws)):
              if isinstance(bw, six.string_types):
                  bw_func = getattr(smnp.bandwidths, "bw_" + bw)
                  bws[i] = bw_func(coord)
          kde = smnp.KDEMultivariate([x, y], "cc", bws)
          x_support = np.linspace(xyranges[0][0], xyranges[0][1], xsize)
          y_support = np.linspace(xyranges[1][0], xyranges[1][1], ysize)
          xx, yy = np.meshgrid(x_support, y_support)
          z = kde.pdf([xx.ravel(), yy.ravel()]).reshape(xx.shape)
          return x_support, y_support, z
docs: Compute a bivariate kde using statsmodels.
          This function is mainly motivated through
          seaborn.distributions._statsmodels_bivariate_kde