semiMixLocalLogLikelihood.hh

Go to the documentation of this file.

#ifndef NPSTAT_SEMIMIXLOCALLOGLIKELIHOOD_HH_
#define NPSTAT_SEMIMIXLOCALLOGLIKELIHOOD_HH_
 
/*!
// \file semiMixLocalLogLikelihood.hh
//
// \brief Local likelihood for fitting semiparametric mixture models
//
// Author: I. Volobouev
//
// August 2022
*/
 
#include <cmath>
#include <vector>
#include <cassert>
#include <climits>
#include <stdexcept>
 
#include "npstat/nm/SimpleFunctors.hh"
#include "npstat/nm/coordAndWeight.hh"
#include "npstat/nm/fcnOrConst.hh"
 
#include "npstat/stat/AbsDistribution1D.hh"
 
namespace npstat {
    // WFcn can be just a constant (convertible to double)
    template <typename Numeric, typename WFcn>
    inline double semiMixLocalLogLikelihood(
        const std::vector<Numeric>& sortedCoords,
        const AbsDistribution1D& signal,
        const double alpha,
        const Functor1<double, double>& bgDensityFcn,
        const WFcn& localizingWeight,
        const double localizingWeightXmin,
        const double localizingWeightXmax)
    {
        if (localizingWeightXmin >= localizingWeightXmax) throw std::invalid_argument(
            "In npstat::semiMixLocalLogLikelihood: invalid weight boundaries");
 
        const double logMin = std::log(DBL_MIN);
 
        long double sum = 0.0L;
        const unsigned long sz = sortedCoords.size();
        double x, dummy;
        for (unsigned long i=0; i<sz; ++i)
        {
            Private::coordAndWeight(sortedCoords[i], &x, &dummy);
 
            // In what follows, the point weights (p.second) are ignored.
            // We assume that those weights were used to solve the Fredholm
            // equation and not imposed by some reason beyond our control.
            if (x >= localizingWeightXmin && x <= localizingWeightXmax)
            {
                const double w = fcnOrConst(localizingWeight, x);
                assert(w >= 0.0);
                if (w > 0.0)
                {
                    const double s = signal.density(x);
                    const double bg = bgDensityFcn(x);
                    const double d = alpha*s + (1.0-alpha)*bg;
                    const double logd = d > DBL_MIN ? std::log(d) : logMin;
                    sum += w*logd;
                }
            }
        }
        return sum;
    }
}
 
#endif // NPSTAT_SEMIMIXLOCALLOGLIKELIHOOD_HH_