Bootstrap.h

Go to the documentation of this file.

// Honeycomb, Copyright (C) 2015 NewGamePlus Inc.  Distributed under the Boost Software License v1.0.
#pragma once

#include "Honey/Math/Random/Random.h"

namespace honey
{


template<class SampleT, sdt Dim = 1, class Real__ = Real>
class Bootstrap
{
public:
    typedef Real__                              Real;
private:
    typedef typename Numeral<Real>::Real_       Real_;
    typedef typename Real_::DoubleType          Double_;
    typedef typename Double_::Real              Double;
    typedef Alge_<Real>                         Alge;
    typedef Interp_<Real>                       Interp;
    typedef Gaussian_<Real>                     Gaussian;

public:
    static const sdt dim                        = Dim;
    typedef Vec<dim,Real>                       Vec;
    typedef vector<const SampleT*>              SampleList;
    typedef function<Vec (const SampleList&)>   Func;


    Bootstrap(const Func& func, RandomGen& gen, const vector<SampleT>& samples, Real alpha = 0.05, szt bootSampleCount = 50) :
        func(func), gen(gen), samples(samples), alpha(alpha), bootSampleCount(bootSampleCount),
        _lower(0), _upper(0), _progress(0), sampleTotal(0), idx(0),
        bootRes(nullptr), jackRes(nullptr), jackMean(0)
    {
        if (samples.size() == 0) return;
        bootRes.set(new Vec[bootSampleCount]);
        bootSamples.resize(samples.size());
        jackRes.set(new Vec[samples.size()]);
    }

    ~Bootstrap() {}


    void calc(Real progressDelta = 1);

    Real progress() const                                                       { return _progress; }
    const Vec& lower() const                                                    { return _lower; }
    const Vec& upper() const                                                    { return _upper; }

    struct meanFunc
    {
        Vec operator()(const SampleList& samples)
        {
            Vec sum = reduce(samples, Vec().fromZero(), [](const Vec& a, const SampleT* e) { return a + *e; });
            return samples.size() > 0 ? sum / (Real)samples.size() : Vec().fromZero();
        }
    };

    struct varianceFunc
    {
        Vec operator()(const SampleList& samples)
        {
            Vec mean = meanFunc()(samples);
            Vec sumDev = reduce(samples, Vec().fromZero(), [&](const Vec& a, const SampleT* e) -> Real { return a + (*e - mean).elemSqr(); });
            return samples.size() > 1 ? sumDev / (Real)(samples.size()-1) : Vec().fromZero();
        }
    };

private:
    //Input
    Func func;
    RandomGen& gen;
    const vector<SampleT>& samples;
    Real alpha;
    szt bootSampleCount;

    //Output
    Vec _lower;
    Vec _upper;

    //Progress
    Real _progress;
    sdt sampleTotal;
    sdt idx;

    //Calc locals
    UniquePtr<Vec[]> bootRes;
    SampleList bootSamples;
    Vec origRes;
    UniquePtr<Vec[]> jackRes;
    Vec jackMean;
};

template<class SampleT, sdt Dim, class Real>
void Bootstrap<SampleT,Dim,Real>::calc(Real progressDelta)
{
    if (_progress == 1)
        return;

    _progress = Alge::min(_progress + progressDelta, 1);
    sdt sampleAcc = _progress*(bootSampleCount + samples.size()) - sampleTotal;
    sampleTotal += sampleAcc;

    if (samples.size() == 0)
        return;

    sdt sampleCount = 0;

    if (sampleTotal - sampleAcc < bootSampleCount)
    {
        //Get function result with bootstrap samples (random sampling with replacement)
        for (; idx < bootSampleCount; ++idx, ++sampleCount)
        {
            if (sampleCount >= sampleAcc)
                return;
            //Build bootstrap samples
            for (szt i = 0; i < samples.size(); ++i)
                bootSamples[i] = &samples[Discrete_d(gen, 0, samples.size()-1).next()];
            //Add bootstrap result
            bootRes[idx] = func(const_cast<const SampleList&>(bootSamples));
        }

        //Get function result with original sample data
        for (szt i = 0; i < samples.size(); ++i)
            bootSamples[i] = &samples[i];
        origRes = func(const_cast<const SampleList&>(bootSamples));

        idx = 0;
    }

    //Get function result with jackknife samples (every sample is omitted once)
    //and calc the jackknife estimator (mean), which is used to estimate bias/variance
    for (; idx < samples.size(); ++idx, ++sampleCount)
    {
        if (sampleCount >= sampleAcc)
            return;
        //Erase one sample
        bootSamples.erase(bootSamples.begin() + idx);
        //Add jackknife result
        jackRes[idx] = func(const_cast<const SampleList&>(bootSamples));
        jackMean += jackRes[idx];
        //Re-insert sample
        bootSamples.insert(bootSamples.begin() + idx, &samples[idx]);
    }
    jackMean /= Real(samples.size());

    vector<Real> ecdf(bootSampleCount);

    for (sdt i = 0; i < Vec::s_size; ++i)
    {
        //Bias, also build empirical cdf (sorted boot results)
        Real sumLess = 0;
        Real sumEq = 0;
        for (szt j = 0; j < bootSampleCount; ++j)
        {
            Real val = bootRes[j][i];
            if (val < origRes[i])
                sumLess++;
            else if (val == origRes[i])
                sumEq++;
            ecdf[j] = val;
        }
        Real z = Gaussian().cdfInv((sumLess + sumEq/2) / bootSampleCount);
        std::sort(ecdf.begin(), ecdf.end());

        //Acceleration
        Real sumDevSqr = 0;
        Real sumDevCube = 0;
        for (szt j = 0; j < samples.size(); ++j)
        {
            Real dev = jackMean[i] - jackRes[j][i];
            Real sqr = Alge::sqr(dev);
            sumDevSqr += sqr;
            sumDevCube += sqr*dev;
        }
        Real acc = (sumDevCube / Alge::pow(sumDevSqr, 1.5)) / 6;

        //Apply bias+acc correction
        Real za = Gaussian().cdfInv(alpha / 2);
        Real idxLower = Gaussian().cdf(z + (z+za) / (1-acc*(z+za))) * (ecdf.size()-1);
        Real idxUpper = Gaussian().cdf(z + (z-za) / (1-acc*(z-za))) * (ecdf.size()-1);
        
        //Empirical cdf, interpolate fractional index
        _lower[i] = idxLower < ecdf.size()-1 ? Interp::linear(Alge::frac(idxLower), ecdf[idxLower], ecdf[idxLower+1]) : ecdf[idxLower];
        _upper[i] = idxUpper < ecdf.size()-1 ? Interp::linear(Alge::frac(idxUpper), ecdf[idxUpper], ecdf[idxUpper+1]) : ecdf[idxUpper];
    }

    //Calculation is complete
    _progress = 1;
}

}