我看过其他关于堆栈交换的例子，但无法拼凑出在我的情况下不起作用的地方，下面是演示问题的精简代码
编辑：我已经更新了FFTWrapper类来进行归一化(通过执行正向和反向序列来校正FFTW应用的nFT因子(
#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;
const double m_sqrt_nFFT;
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_sqrt_nFFT(std::sqrt(m_nFFT))
{
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]/m_sqrt_nFFT;
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
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);
std::vector<double> out(m_nFFT);
for(i=0; i<m_nFFT; ++i)
out[i]=m_real[i]/m_sqrt_nFFT;
return out;
}
};
// 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;
}
// iChange[nC-1][0]=0.0; //some fully complex examples have this but it doesn't appear to make a difference
// iChange[nC-1][1]=0.0;
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));
std::vector<double> cdf2(fftw.inverse(f));
inPlaceFFTDerivative(f, nFFT, dx);
std::vector<double> pdf(fftw.inverse(f));
for(i=0; i<nSegEdgesPerInstance; ++i)
std::cerr << pEdges[i] << " " << cdf[i] << " vs " << cdf2[i] << ": " << pdfShouldBe[i] << " vs " << pdf[i] << std::endl;
return 0;
}
这是的输出
0对1.81299e-16:0对-11.1996
0.0625 1.77057e-13对1.77401e-13:3.17909e-11对4.93837
0.125 4.16583e-10对4.16583e-10:3.7231e-08对-3.11885
0.1875 3.8221e-08对3.8221e-08:2.26554e-06对2.27159
0.25 9.26922e-07对9.26922e-07:4.99623e-05对-1.78611
0.3125 1.08197e-05对1.08197e-05:0.000379853对1.47352
0.375 7.93519e-05 vs 7.93519e-05:0.00230244 vs-1.25323
0.4375 0.0004213 vs 0.0004213：0.0103742 vs 1.10674
0.5 0.00176101 vs 0.00176101：0.0374823 vs-0.938405
0.5625 0.00611262 vs 0.00611262：0.113884 vs 0.996605
0.625 0.0182524 vs 0.0182524：0.3000018 vs-0.510037
0.6875 0.0480048 vs 0.0480048：0.698306 vs 1.4521
0.75 0.112878 vs 0.112878：1.44857 vs 0.737019
0.8125 0.238987 vs 0.238987：2.66815 vs 3.3495
0.8750.454276 vs 0.454276:4.24137 vs 3.5857
0.9375 0.757705 vs 0.757705:5.18051 vs 5.74862
1 1 vs 1：0 vs 1.22598
我的问题是：我到底把这个FFTW3衍生产品做错了什么？
潜在的问题是fft假设函数是周期性
PDF开始(在-inf或更晚(，结束(在+inf或更早(在yvalue=0，因此您可以将它们包裹起来/视为周期性的。
CDF以yvalue=0开始，以yvalue=1结束，因此不能将其视为周期性
正因为如此，这种尝试注定会失败，双环卷积可以使用CDF进行，但那时你无法获得漂亮的FFT导数。