diff --git a/PeakOTron.py b/PeakOTron.py index 22c6145fcef2aced8aafeb1c45e9e2cf19c99840..0731dcd172bc79be238616b3a7137fa7d3f692b1 100644 --- a/PeakOTron.py +++ b/PeakOTron.py @@ -49,6 +49,9 @@ class FakeFuncCode: class BandWidthOptimiser: + #Hobbema, Hermans, and Van den Broeck (1971) and by Duin (1976) + #leave-one-out cross-validation approach + #https://cran.r-project.org/web/packages/kedd/vignettes/kedd.pdf def _Dummy(self, _, bw): return bw @@ -56,9 +59,12 @@ class BandWidthOptimiser: def _KDE(self, bw): try: x_kde, y_kde = FFTKDE(kernel = self.kernel, bw=bw).fit(self.data).evaluate(self.n_kde_samples) - return np.log(y_kde) + + loss = np.log((self.N-1)*bw) - np.sum(np.log(y_kde)) + + return loss except: - return np.empty(self.n_kde_samples) + return np.nan def _PPF(self, data): """ Compute ECDF """ @@ -77,6 +83,7 @@ class BandWidthOptimiser: self.data = data self.PPF = self._PPF(self.data) self.data = self.data[(self.data<self.PPF(alpha))] + self.N = len(data) self.kernel=kernel self.n_kde_samples=n_kde_samples self.last_arg = None @@ -86,13 +93,10 @@ class BandWidthOptimiser: def __call__(self, *arg): self.last_arg = arg - log_y_hat = self._KDE(*arg) - loss = -np.sum(log_y_hat) - + loss = self._KDE(*arg) return loss - class BinnedLH: @@ -250,6 +254,7 @@ class PeakOTron: "omega"] self._eps = np.finfo(np.float64).eps * 10 + self._eps_kde = 1e-6 self._FWHM2Sigma = 1/(2*np.sqrt(2*np.log(2))) self._plot_figsize= (10,10) @@ -482,8 +487,8 @@ class PeakOTron: try: BWO= BandWidthOptimiser(data, kernel=kernel) bw_scott = scott_bin_width(data) - minuit_kde = Minuit(BWO, bw=bw_scott) - minuit_kde.limits["bw"]=(self._eps, bw_scott) + minuit_kde = Minuit(BWO, bw=bw_scott/2) + minuit_kde.limits["bw"]=(self._eps_kde, 3*bw_scott) minuit_kde.strategy=2 minuit_kde.migrad(ncall= self._n_call_minuit, iterate =self._n_iterations_minuit) @@ -2087,7 +2092,7 @@ class PeakOTron: print("Performing fit initialisation...") self.Init() self.InitFit(data, **kwargs_fit) - + if(not self._failed): print("Fitting...")