概率密度的正交例程

Question

我想整合概率密度函数,(-\infty, a]因为cdf不是封闭形式.但我不确定如何在C++中这样做.

这个任务在Mathematica中非常简单; 我需要做的就是定义函数,

f[x_, lambda_, alpha_, beta_, mu_] := 
   Module[{gamma}, 
     gamma = Sqrt[alpha^2 - beta^2]; 
     (gamma^(2*lambda)/((2*alpha)^(lambda - 1/2)*Sqrt[Pi]*Gamma[lambda]))*
      Abs[x - mu]^(lambda - 1/2)*
      BesselK[lambda - 1/2, alpha Abs[x - mu]] E^(beta (x - mu))
   ];

然后调用NIntegrate例程以数字方式集成它.

F[x_, lambda_, alpha_, beta_, mu_] := 
    NIntegrate[f[t, lambda, alpha, beta, mu], {t, -\[Infinity], x}] 

现在我想在C++中实现同样的功能.我使用gsl_integration_qagilgsl数字库中的例程.它旨在将功能集成在半无限区间(-\infty, a],这正是我想要的.但不幸的是,我无法让它发挥作用.

这是C++中的密度函数,

density(double x)
{
using namespace boost::math;

if(x == _mu)
    return std::numeric_limits<double>::infinity();

    return pow(_gamma, 2*_lambda)/(pow(2*_alpha, _lambda-0.5)*sqrt(_pi)*tgamma(_lambda))* pow(abs(x-_mu), _lambda - 0.5) * cyl_bessel_k(_lambda-0.5, _alpha*abs(x - _mu)) * exp(_beta*(x - _mu));

}  

然后我尝试通过调用gsl例程来集成以获取cdf.

cdf(double x)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

    double result, error;      
    gsl_function F;
    F.function = &density;

    double epsabs = 0;
    double epsrel = 1e-12;

    gsl_integration_qagil (&F, x, epsabs, epsrel, 1000, w, &result, &error);

    printf("result          = % .18f\n", result);
    printf ("estimated error = % .18f\n", error);
    printf ("intervals =  %d\n", w->size);

    gsl_integration_workspace_free (w);

    return result;

}

但是会gsl_integration_qagil返回错误number of iterations was insufficient.

 double mu = 0.0f;
 double lambda = 3.0f;
 double alpha = 265.0f;
 double beta = -5.0f;

 cout << cdf(0.01) << endl;

如果我增加工作空间的大小,那么bessel函数将不会评估.

我想知道是否有人可以让我对我的问题有所了解.的呼叫到相应的数学函数F以上x = 0.01的回报0.904384.

可能是密度集中在一个非常小的区间(即[-0.05, 0.05]密度几乎在外面,0下面给出了一个图).如果是这样的话可以做些什么.谢谢阅读.

Answer 1

回复：积分到+/-无穷大：

我将使用 Mathematica 来找到 |x - μ| 的经验界限 >> K，其中 K 表示平均值周围的“宽度”，K 是 alpha、beta 和 lambda 的函数 - 例如 F 小于并约等于 a(x-μ) ^-2或 ae ^{- b(x-μ) 2}或其他。这些函数具有已知的无穷大积分，您可以根据经验对其进行评估。然后，您可以对 K 进行数值积分，并使用有界近似从 K 到无穷大。

计算出 K 可能有点棘手；我对贝塞尔函数不太熟悉，所以我无法为您提供太多帮助。

一般来说，我发现对于并不明显的数值计算，最好的方法是在进行数值评估之前尽可能多地进行分析数学。（有点像自动对焦相机 - 将其靠近您想要的位置，然后让相机完成其余的工作。）