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``.

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()