sampleDistro1DWithWeight.hh

Go to the documentation of this file.

#ifndef NPSTAT_SAMPLEDISTRO1DWITHWEIGHT_HH_
#define NPSTAT_SAMPLEDISTRO1DWITHWEIGHT_HH_
 
/*!
// \file sampleDistro1DWithWeight.hh
//
// \brief Acceptance-rejection sampling from a density with compensatory
//        weights. Does not change the average but affects the variance.
//
// Author: I. Volobouev
//
// June 2022
*/
 
#include <cmath>
#include <vector>
#include <utility>
#include <cassert>
#include <stdexcept>
#include <climits>
 
#include "npstat/nm/Triple.hh"
#include "npstat/nm/fcnOrConst.hh"
 
#include "npstat/rng/AbsRandomGenerator.hh"
#include "npstat/stat/AbsDistribution1D.hh"
 
namespace npstat {
    enum SampleSizeInterpretation
    {
        SIZE_IS_N_TRIES = 0,
        SIZE_IS_N_GENERATED,
        SIZE_IS_KISHS
    };
 
    /**
    // AcceptanceFunction is a functor which takes a double as
    // an argument and returns a double on the [0, 1] interval.
    //
    // The elements of the returned triple are:
    //
    //  first  -- sum of weights
    //
    //  second -- the Kish's effective sample size of the generated
    //            sample
    //
    //  third  -- the efficiency of the generator (the ratio between
    //            the number of points generated and the number of
    //            points attempted)
    //
    // In each generated point (the vector "sample" is filled with them),
    // the first element of the pair is the point coordinate and the
    // second element is the point weight.
    */
    template <class AcceptanceFunction>
    inline npstat::Triple<double,double,double> sampleDistro1DWithWeight(
        const AbsDistribution1D& distro,
        const AcceptanceFunction& fcn,
        AbsRandomGenerator& rng,
        const double sampleSize,
        const SampleSizeInterpretation whichSize,
        std::vector<std::pair<double,double> >* sample)
    {
        if (sampleSize < 0.0) throw std::invalid_argument(
            "In npstat::sampleDistro1DWithWeight: "
            "sample size can not be negative");
        if (sample)
            sample->clear();
        if (sampleSize > 0.0)
        {
            assert(sample);
 
            long double kishLimit = LDBL_MAX;
            unsigned long nTriesLimit = ULONG_MAX;
            unsigned long nGeneratedLimit = ULONG_MAX;
            double lastTryProb = 1.0;
            double lastGenProb = 1.0;
 
            switch (whichSize)
            {
            case SIZE_IS_N_TRIES:
                nTriesLimit = std::floor(sampleSize);
                if (static_cast<double>(nTriesLimit) != sampleSize)
                {
                    lastTryProb = sampleSize - nTriesLimit;
                    ++nTriesLimit;
                }
                break;
 
            case SIZE_IS_N_GENERATED:
                nGeneratedLimit = std::floor(sampleSize);
                if (static_cast<double>(nGeneratedLimit) != sampleSize)
                {
                    lastGenProb = sampleSize - nGeneratedLimit;
                    ++nGeneratedLimit;
                }
                break;
 
            case SIZE_IS_KISHS:
                kishLimit = sampleSize;
                break;
 
            default:
                assert(!"In npstat::sampleDistro1DWithWeight: "
                       "unhandled switch case. This is a bug. Please report.");
            }
 
            unsigned long nTries = 0, nAccepted = 0;
            long double effSize = 0.0L, sumW = 0.0L, sumWSq = 0.0L;
            double tmp;
 
            while (effSize < kishLimit &&
                   nAccepted < nGeneratedLimit &&
                   nTries < nTriesLimit)
            {
                if (nTries + 1U == nTriesLimit)
                    if (rng() > lastTryProb)
                        break;
 
                distro.random(rng, &tmp);
                const double sf = fcnOrConst(fcn, tmp);
                assert(sf >= 0.0);
                assert(sf <= 1.0);
                if (sf > 0.0)
                {
                    const double r = rng();
                    if (r <= sf)
                    {
                        ++nAccepted;
                        bool accepted = true;
                        if (nAccepted == nGeneratedLimit)
                            if (rng() > lastGenProb)
                                accepted = false;
                        if (accepted)
                        {
                            const double w = 1.0/sf;
                            sample->push_back(std::pair<double,double>(tmp, w));
                            sumW += w;
                            sumWSq += w*w;
                            effSize = sumW/sumWSq*sumW;
                        }
                    }
                }
 
                ++nTries;
            }
 
            double eff = 0.0;
            if (nTries)
                eff = static_cast<double>(nAccepted)/nTries;
            return npstat::Triple<double,double,double>(sumW, effSize, eff);
        }
        else
            return npstat::Triple<double,double,double>(0.0, 0.0, 0.0);
    }
}
 
#endif // NPSTAT_SAMPLEDISTRO1DWITHWEIGHT_HH_