docs: def _make_kde(self, use_sklearn=False, bandwidth=None, rtol=1e-6,
                    sig_clip=50, no_sig_clip=False, cov_all=True,
                    **kwargs):
          """Creates KDE objects for 3-d shape parameter distribution
  
          KDE represents likelihood as function of trapezoidal
          shape parameters (log(delta), T, T/tau).
  
          Uses :class:`scipy.stats.gaussian_kde`` KDE by default;
          Scikit-learn KDE implementation tested a bit, but not
          fully implemented.
  
          :param use_sklearn:
              Whether to use scikit-learn implementation of KDE.
              Not yet fully implemented, so this should stay ``False``.
  
          :param bandwidth, rtol:
              Parameters for sklearn KDE.
  
          :param **kwargs:
              Additional keyword arguments passed to
              :class:`scipy.stats.gaussian_kde``.
  
          """
  
          try:
              #define points that are ok to use
              first_ok = ((self.stars['slope'] > 0) &
                          (self.stars['duration'] > 0) &
                          (self.stars['duration'] < self.period) &
                          (self.depth > 0))
          except KeyError:
              logging.warning('Must do trapezoid fits before making KDE.')
              return
  
          self.empty = False
          if first_ok.sum() < 4:
              logging.warning('Empty population ({}): < 4 valid systems! Cannot calculate lhood.'.format(self.model))
              self.is_empty = True #will cause is_ruled_out to be true as well.
              return
              #raise EmptyPopulationError('< 4 valid systems in population')
  
          logdeps = np.log10(self.depth)
          durs = self.stars['duration']
          slopes = self.stars['slope']
  
          #Now sigma-clip those points that passed first cuts
          ok = np.ones(len(logdeps), dtype=bool)
          for x in [logdeps, durs, slopes]:
              med = np.median(x[first_ok])
              mad = np.median(np.absolute(x[first_ok] - med))
              ok &= np.absolute(x - med) / mad < sig_clip
  
          second_ok = first_ok & ok
  
          # Before making KDE for real, first calculate
          #  covariance and inv_cov of uncut data, to use
          #  when it's cut, too.
  
          points = np.array([durs[second_ok],
                             logdeps[second_ok],
                             slopes[second_ok]])
          kde = gaussian_kde(np.vstack(points)) #backward compatibility?
          cov_all = kde._data_covariance
          icov_all = kde._data_inv_cov
          factor = kde.factor
  
          # OK, now cut the data for constraints & proceed
  
          ok = second_ok & self.distok
  
          logdeps = logdeps[ok]
          durs = durs[ok]
          slopes = slopes[ok]
  
          if ok.sum() < 4 and not self.empty:
              logging.warning('Empty population ({}): < 4 valid systems! Cannot calculate lhood.'.format(self.model))
              self.is_empty = True
              return
              #raise EmptyPopulationError('< 4 valid systems in population')
  
  
          if use_sklearn:
              self.sklearn_kde = True
              logdeps_normed = (logdeps - logdeps.mean())/logdeps.std()
              durs_normed = (durs - durs.mean())/durs.std()
              slopes_normed = (slopes - slopes.mean())/slopes.std()
  
              #TODO: use sklearn preprocessing to replace below
              self.mean_logdepth = logdeps.mean()
              self.std_logdepth = logdeps.std()
              self.mean_dur = durs.mean()
              self.std_dur = durs.std()
              self.mean_slope = slopes.mean()
              self.std_slope = slopes.std()
  
              points = np.array([logdeps_normed, durs_normed, slopes_normed])
  
              #find best bandwidth.  For some reason this doesn't work?
              if bandwidth is None:
                  grid = GridSearchCV(KernelDensity(rtol=rtol),
                                      {'bandwidth':np.linspace(0.05,1,50)})
                  grid.fit(points)
                  self._best_bandwidth = grid.best_params_
                  self.kde = grid.best_estimator_
              else:
                  self.kde = KernelDensity(rtol=rtol, bandwidth=bandwidth).fit(points)
          else:
              self.sklearn_kde = False
              points = np.array([durs, logdeps, slopes])
              self.kde = gaussian_kde(np.vstack(points), **kwargs) #backward compatibility?
  
              # Reset covariance based on uncut data
              self.kde._data_covariance = cov_all
              self.kde._data_inv_cov = icov_all
              self.kde._compute_covariance()