how do I take a 1st derivative of a CDF with FFTW3 r2c

  c++, cdf, derivative, fftw, pdf

I’ve looked at other examples on stack exchange and haven’t been able to piece together what doesn’t work in my case, here is the stripped down code demonstrating the issue

#include <cmath>
#include <fftw3.h>
#include <vector>
#include <boost/math/special_functions/beta.hpp>

//! Danger: if you move (i.e. copy and delete original) you'll get a memory error, so use a shared_ptr around FFTWrapper and you'll be safe
class FFTWrapper
{
public:
    const unsigned m_nFFT;
    const unsigned m_nC;
protected:
    fftw_plan m_forward;
    fftw_plan m_inverse;
    mutable std::vector<double> m_real;
    mutable fftw_complex* m_complex;
public:

    FFTWrapper(const unsigned nFFT) : m_nFFT(nFFT), m_nC((m_nFFT / 2 + 1))
    {
        m_real.resize(nFFT);
        m_complex = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * m_nC);
        m_forward = fftw_plan_dft_r2c_1d(m_nFFT, &(m_real[0]), m_complex, FFTW_ESTIMATE);
        m_inverse = fftw_plan_dft_c2r_1d(m_nFFT, m_complex, &(m_real[0]), FFTW_ESTIMATE);
    };

    ~FFTWrapper()
    {
        fftw_destroy_plan(m_forward);
        fftw_destroy_plan(m_inverse);
        fftw_free(m_complex);
    };

    inline unsigned nFftReal() const { return m_nFFT; };
    inline unsigned nFftComplex() const { return m_nC; };

    std::vector<std::vector<double>> forward(const std::vector<double>& in)
    {
        unsigned i(0);
        assert(in.size()==m_nFFT);
        for(; i<m_nFFT; ++i)
            m_real[i]=in[i];
        fftw_execute(m_forward);
        std::vector<std::vector<double>> out(m_nC,std::vector<double>(2));
        for(i=0; i<m_nC; ++i)
        {
            out[i][0]=m_complex[i][0];
            out[i][1]=m_complex[i][1];
        }
        return out;
    };

    //! copy the output because it will change when eeiterh forward() or inverse() is called
    const std::vector<double>& inverse(const std::vector<std::vector<double>>& in)
    {
        unsigned i(0);
        assert((in[0].size()==2)&&(in.size()==m_nC));
        for(; i<m_nC; ++i)
        {
            m_complex[i][0]=in[i][0];
            m_complex[i][1]=in[i][1];
        }
        fftw_execute(m_inverse);
        return m_real;
    }
};

// FFTW3 documentation says it works best when nFFT is a product of small prime numbers
// this "rounds up" nFFT to a product of mostly 2's (for more accurate division), 
// optionally one or two 3's, and optionally one 5
unsigned niceNFFT(const unsigned number)
{
    const double nl2(std::max(2.0,std::ceil(std::log(number)/std::log(2.0))));
    if(nl2<=6.0)
        return std::pow(2.0,nl2);
    if(number<=96)
        return 96;
    if(number<=128)
        return 128;
    const unsigned f(std::pow(2.0,nl2-8.0));
    if(number<=f*144)
        return f*144;
    if(number<=f*160)
        return f*160;
    if(number<=f*192)
        return f*192;
    if(number<=f*240)
        return f*240;
    if(number<=f*256)
        return f*256;
    assert(false);
}

inline std::vector<std::vector<double>>& inPlaceFFTDerivative(std::vector<std::vector<double>>& iChange, const double nFFT, const double dx)
{
    const unsigned nC(iChange.size());
    if(nC==0)
        return iChange;
    assert(nC==std::floor(nFFT/2.0)+1);
    assert(iChange[0].size()==2);
    double dblTemp, kappa; 
    for(unsigned j=0; j<nC; ++j)
    {
        kappa=2.0*M_PI*static_cast<double>(j)/(nFFT*dx);
        dblTemp=-iChange[j][1]*kappa;
        iChange[j][1]=iChange[j][0]*kappa;
        iChange[j][0]=dblTemp;
    }
    return iChange;
}

//assuming binomial draws -> nSucc & nFail which we want to use estimate the probability for this instance
//the beta distribution is the Bayesian conjugate prior for the binomial distribution
void posteriorBetaCdf(std::vector<double>& cdf, 
    const std::vector<double>& pEdges, //!< these must be evenly spaced between 0 and 1 (inclusive)
    const double alpha, const double beta)
{
    const unsigned nSegEdges(pEdges.size());
    assert(nSegEdges<=cdf.size());
    cdf[0]=0.0; //question should cdf[0] always be zero? ibeta(0,n,0.0) says 1.0, matlab's betainc(0.0,0,n) says 0.0
    for(unsigned i=1; i<nSegEdges-1; ++i)
        cdf[i]=boost::math::ibeta(alpha,beta,pEdges[i]); //this function is SLOW!!!!!!!
    for(unsigned i=nSegEdges-1; i<cdf.size(); ++i)
        cdf[i]=1.0;
}

//this stripped down (for StackOverflow) version of function, doesn't handle alpha<1 or beta<1 cases
void posteriorBetaPdf(std::vector<double>& pdf, 
    const std::vector<double>& pEdges, //!< these must be evenly spaced between 0 and 1 (inclusive)
    const double alpha, const double beta)
{
    const unsigned nSegEdges(pEdges.size());
    assert(nSegEdges<=pdf.size());
    double am1(alpha-1.0), bm1(beta-1.0), denom(boost::math::beta(alpha,beta));
    for(unsigned i=0; i<nSegEdges; ++i)
        pdf[i]=std::pow(pEdges[i],am1)*std::pow(1.0-pEdges[i],bm1)/denom;
    for(unsigned i=nSegEdges; i<pdf.size(); ++i)
        pdf[i]=0.0;
}


//intent is to EVENTUALLY convolve cdfs (rather than pdfs for numerical advantages) using FFTs AND
//take FFT derivative each time to always have the FFT of a CDF, but I'm getting stuck on getting 
//the FFT derivative working
int main()
{
    const unsigned nSample(15), nSucc(13), nFail(2);
    const unsigned nSegsPerInstance(16); //be kind to computer by rounding nSample up to a power of 2
    const unsigned nSegEdgesPerInstance(nSegsPerInstance+1);
    const unsigned nInstancesToSum(5);  //would be larger in real application
    const unsigned nSegsTotal(nSegsPerInstance*nInstancesToSum);
    const unsigned nSegEdgesTotal(nSegsTotal+1);
    const unsigned nFFT(niceNFFT(nSegEdgesTotal));
    double dx(1.0/nSegsPerInstance); //only correct before the first convolution 
    std::vector<double> cdf(nFFT,1.0), pdfShouldBe(nFFT,0.0), pEdges(nSegsPerInstance+1);
    double deBias(static_cast<double>(nSample-2)/static_cast<double>(nSample)); //this makes beta posterior mean and stddev equal the sample mean and stddev (of an indicator function)
    double alpha(nSucc*deBias), beta(nFail*deBias);

    unsigned i(0);
    for(; i<nSegsPerInstance; ++i)
        pEdges[i]=i*dx;
    pEdges[nSegsPerInstance]=1.0;

    posteriorBetaCdf(cdf,pEdges,alpha,beta);
    posteriorBetaPdf(pdfShouldBe,pEdges,alpha,beta);

    FFTWrapper fftw(nFFT);

    std::vector<std::vector<double>> f(fftw.forward(cdf));
    inPlaceFFTDerivative(f, nFFT, dx);
    std::vector<double> pdf(fftw.inverse(f));

    for(i=0; i<nSegEdgesPerInstance; ++i)
        std::cerr << pEdges[i] << " " << cdf[i] << ": " << pdfShouldBe[i] << " vs " << pdf[i] << std::endl;

    return 0;
}

and here’s the output

pEdges cdf: pdfShouldBe vs pdf
0 0: 0 vs -1067.48
0.0625 1.77057e-13: 3.17909e-11 vs 474.083
0.125 4.16583e-10: 3.7231e-08 vs -299.41
0.1875 3.8221e-08: 2.26554e-06 vs 218.073
0.25 9.26922e-07: 4.09623e-05 vs -171.467
0.3125 1.08197e-05: 0.000379853 vs 141.458
0.375 7.93519e-05: 0.00230244 vs -120.31
0.4375 0.0004213: 0.0103742 vs 106.247
0.5 0.00176101: 0.0374823 vs -90.0869
0.5625 0.00611262: 0.113884 vs 95.6741
0.625 0.0182524: 0.300018 vs -48.9635
0.6875 0.0480048: 0.698306 vs 139.402
0.75 0.112878: 1.44857 vs 70.7538
0.8125 0.238987: 2.66815 vs 321.552
0.875 0.454276: 4.24137 vs 344.227
0.9375 0.757705: 5.18051 vs 551.867
1 1: 0 vs 117.694

My question is: what the heck am I doing wrong with this FFTW3 derivative?

Source: Windows Questions C++

LEAVE A COMMENT