amiseOptimalBandwidth.hh

Go to the documentation of this file.

#ifndef NPSTAT_AMISEOPTIMALBANDWIDTH_HH_
#define NPSTAT_AMISEOPTIMALBANDWIDTH_HH_
 
/*!
// \file amiseOptimalBandwidth.hh
//
// \brief Optimal AMISE bandwidth for KDE with high order kernels
//
// AMISE stands for Asymptotic Mean Integrated Squared Error.
// The formulae used in this code come from the paper "Bandwidth
// Selection in Kernel Density Estimation: A Review" by B.A. Turlach.
//
// Author: I. Volobouev
//
// July 2010
*/
 
namespace npstat {
    /**
    // Estimate optimal bandwidth for filters generated by the Gaussian
    // distribution used as the weight function. The arguments are as follows:
    //
    // filterDegree  -- Degree of the filter (in the LOrPE sense). This
    //                  degree is less by 2 than the order of the kernel
    //                  in the Turlach's paper. For example, filter of
    //                  degree 0 corresponds to kernel of order 2 (original
    //                  Gaussian), filter of degree 2 corresponds to kernel
    //                  of order 4, etc. Here, the degree must be even.
    //                  Maximum possible filter degree is 42.
    //
    // npoints       -- Number of points in the data sample.
    //
    // fvalues       -- Array of scanned values of the reference density.
    //                  Note that the contents of this array are NOT preserved.
    //
    // arrLen        -- Number of elements in the array "fvalues".
    //
    // h             -- Step with which the scan of the reference density
    //                  was performed.
    //
    // expectedAmise -- If this argument is provided, it will be filled
    //                  with the expected AMISE value.
    //
    // The formulae used involve numerical evaluations of the derivatives
    // of the scanned reference density. For a given filterDegree, the
    // code needs to calculate the derivative of order filterDegree + 2.
    // Make sure the scan step is chosen appropriately for this kind of
    // calculation (see, for example, the "Numerical Recipes" book in order
    // to understand the issues involved).
    */
    template<typename Real>
    double amiseOptimalBwGauss(unsigned filterDegree, double npoints,
                               Real* fvalues, unsigned long arrLen,
                               Real h, double* expectedAmise = 0);
 
    /**
    // Estimate optimal bandwidth for filters generated by the symmetric
    // beta distribution used as the weight function. The "power" argument
    // is the power of the symmetric beta, limited here to unsigned integers.
    // Maximum possible power value is 10. The meaning of all other arguments
    // is the same as in the function npstat::amiseOptimalBwGauss.
    */
    template<typename Real>
    double amiseOptimalBwSymbeta(unsigned power, unsigned filterDegree,
                                 double npoints,
                                 Real* fvalues, unsigned long arrLen,
                                 Real h, double* expectedAmise = 0);
 
    /**
    // Plug-in bandwidth for filters generated by the Gaussian distribution
    // used as the weight function
    */
    double amisePluginBwGauss(unsigned filterDegree, double npoints,
                              double sampleSigma, double* expectedAmise = 0);
 
    /**
    // Approximate version of "amisePluginBwGauss" for use with non-integer
    // filter degrees
    */
    double approxAmisePluginBwGauss(double filterDegree, double npoints,
                                    double sampleSigma);
 
    /**
    // Plug-in bandwidth for filters generated by the symmetric beta
    // distribution used as the weight function. The "power" argument
    // is the power of the symmetric beta, limited here to unsigned integers.
    // Maximum possible power value is 10.
    */
    double amisePluginBwSymbeta(unsigned power, unsigned filterDegree,
                                double npoints, double sampleSigma,
                                double* expectedAmise = 0);
 
    /**
    // Ratio of the symmetric beta kernel AMISE bandwidth to the Gaussian
    // kernel AMISE bandwidth (as in the concept of "canonical bandwidth").
    //
    // 1.0 will be returned in case the "power" argument is negative.
    */
    double symbetaBandwidthRatio(int power, unsigned filterDegree);
 
    /**
    // Approximate continuous version of "symbetaBandwidthRatio" for use with
    // non-integer filter degrees (uses a polynomial fit to the original).
    // Can be used in combination with "approxAmisePluginBwGauss" to derive
    // a continuous version of symmetric beta bandwidth.
    */
    double approxSymbetaBandwidthRatio(int power, double filterDegree);
 
    /** AMISE-optimal filter degree for the Gaussian kernel */
    unsigned amisePluginDegreeGauss(double npoints);
 
    /** AMISE-optimal filter degree for symmetric beta kernels */
    unsigned amisePluginDegreeSymbeta(unsigned power, double npoints);
 
    /** Maximum filter degree supported by AMISE calculations */
    unsigned maxFilterDegreeSupported();
 
    /** 
    // Integral of kernel squared for Gaussian (power < 0) and
    // symmetric Beta kernels
    */
    double integralOfSymmetricBetaSquared(int power);
 
    /** 
    // Integral of kernel squared for Gaussian (power < 0) and
    // symmetric Beta kernels within some definite limits
    */
    double integralOfSymmetricBetaSquared(int power, double a, double b);
}
 
#ifdef SWIG
 
#include <cassert>
#include <vector>
#include <utility>
 
namespace npstat {
    inline std::pair<double, double>
    amiseOptimalBwGauss2(unsigned filterDegree, double npoints,
                         const double* fvalues, unsigned long arrLen,
                         double h)
    {
        assert(fvalues);
        assert(arrLen);
 
        std::vector<double> tmp(fvalues, fvalues+arrLen);
        double expectedAmise;
        const double bw = amiseOptimalBwGauss(filterDegree, npoints, &tmp[0],
                                              arrLen, h, &expectedAmise);
        return std::make_pair(bw, expectedAmise);
    }
 
    inline std::pair<double, double>
    amiseOptimalBwSymbeta2(unsigned power, unsigned filterDegree,
                           double npoints,
                           const double* fvalues, unsigned long arrLen,
                           double h)
    {
        assert(fvalues);
        assert(arrLen);
 
        std::vector<double> tmp(fvalues, fvalues+arrLen);
        double expectedAmise;
        const double bw = amiseOptimalBwSymbeta(power, filterDegree,
                                                npoints, &tmp[0],
                                                arrLen, h, &expectedAmise);
        return std::make_pair(bw, expectedAmise);
    }
}
#endif // SWIG
 
#include "npstat/stat/amiseOptimalBandwidth.icc"
 
#endif // NPSTAT_AMISEOPTIMALBANDWIDTH_HH_