Bootstrapping.py

from HelperFunctions import FakeFuncCode, SelectRangeNumba, EmpiricalPPF
from KDEpy import FFTKDE
import numpy as np
from scipy.interpolate import interp1d
from iminuit import Minuit
from astropy.stats import bootstrap, scott_bin_width
import scipy.special as sc
import logging

            
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
    
    
    def _KDE(self, bw):

        x_kde, y_kde = FFTKDE(kernel = self.kernel, bw=bw).fit(self.data).evaluate(self.n_kde_samples)        
        loss = -np.mean(np.log(y_kde))
        return loss


        
    def __init__(self, data, kernel="gaussian", alpha = 0.95, n_kde_samples=2**13):
        if(alpha<=0.5):
            raise ValueError("alpha must be above 0.5")
        
        self.data = data
        PPF = EmpiricalPPF(self.data)
        x_low, x_hi = PPF(0.5-alpha/2), PPF(0.5+alpha/2)
        self.data = SelectRangeNumba(self.data, x_low, x_hi)
        self.N = len(data)
        self.kernel=kernel
        self.n_kde_samples=n_kde_samples
        self.last_arg = None
        self.func_code = FakeFuncCode(self._Dummy, dock=True)
        self.eps = np.finfo(np.float64).eps * 10
        self.eps_inv = 1/self.eps
        

    def __call__(self, *arg):
        self.last_arg = arg
        loss = self._KDE(*arg)
        return loss
    
    
def Bootstrap(data, statistic, n_bootstrap, alpha=0.95, ):
    if not (0 < alpha < 1):
        raise ValueError("confidence level must be in (0, 1)")


    if len(data) < 1:
        raise ValueError("data must contain at least one measurement.")


    boot_stat = bootstrap(data, n_bootstrap, bootfunc=statistic)

    stat_data = statistic(data)
    mean_stat = np.mean(boot_stat)
    est_stat = 2*stat_data - mean_stat
    std_err = np.std(boot_stat)
    z_score = np.sqrt(2.0)*sc.erfinv(alpha)
    conf_interval = est_stat + z_score*np.array((-std_err, std_err))


    return est_stat, std_err, conf_interval
    
    
    
    
def BootstrapKDE(data,
                 n_bootstrap,
                 n_call=1000,
                 n_iterations=10,
                 kernel = "gaussian",
                 n_kde_samples=2**14,
                 alpha=0.95,
                 bw_limits=(1e-6, None),
                 limits=None,

                ):

    if not (0 < alpha < 1):
        raise ValueError("Bootstrap confidence level, alpha, must be in (0, 1).")

    if len(data) <= 3:
        raise ValueError("Bootstrap data must contain at least three measurements.")

    try:
        logging.debug("Optimising bandwidth...")
        BWO= BandWidthOptimiser(data, kernel=kernel)
        bw_scott = scott_bin_width(data)
        minuit_kde = Minuit(BWO, bw=bw_scott/2)
        minuit_kde.limits["bw"]=bw_limits
        minuit_kde.strategy=2
        minuit_kde.migrad(ncall=n_call,
                   iterate=n_iterations)

        bw = minuit_kde.values["bw"]
        logging.debug("Bandwidth optimised to {:3.3E} units.".format(bw))
    except:
        logging.warning("Bandwidth optimisation failed. Defaulting to Improved Sheather Jones...")
        bw = "ISJ"


    kde = FFTKDE(kernel = kernel, bw=bw).fit(data)

    x_kde, y_kde_orig = kde.evaluate(n_kde_samples)


    boot_data = bootstrap(data, n_bootstrap)
    y_kde_bs = np.vstack([FFTKDE(kernel = kernel, bw=bw).fit(_data).evaluate(x_kde)
                               for _data in boot_data])


    y_kde_mean = np.mean(y_kde_bs, axis=0)
    y_kde = 2*y_kde_orig - y_kde_mean
    y_kde_err = np.std(y_kde_bs, axis=0)
    z_score = np.sqrt(2.0)*sc.erfinv(alpha)
    y_kde_conf_interval = y_kde + z_score*np.array((-y_kde_err, y_kde_err))

    if(limits is not None):
        cond_inrange =(x_kde>np.min(limits)) & (x_kde<np.max(limits))
        x_kde = x_kde[cond_inrange]
        y_kde = y_kde[cond_inrange]
        y_kde_err = y_kde_err[cond_inrange]

    return x_kde, y_kde, y_kde_err, y_kde_conf_interval


def Bootstrap(data, statistic, n_bootstrap, alpha=0.95):
    if not (0 < alpha < 1):
        raise ValueError("confidence level must be in (0, 1)")


    if len(data) < 1:
        raise ValueError("data must contain at least one measurement.")


    boot_stat = bootstrap(data, n_bootstrap, bootfunc=statistic)

    stat_data = statistic(data)
    mean_stat = np.mean(boot_stat)
    est_stat = 2*stat_data - mean_stat
    std_err = np.std(boot_stat)
    z_score = np.sqrt(2.0)*sc.erfinv(alpha)
    conf_interval = est_stat + z_score*np.array((-std_err, std_err))


    return est_stat, std_err, conf_interval