c++ - 概率密度的正交例程

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 数值库中的例程 gsl_integration_qagil。它旨在集成半无限区间 (-\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 返回错误，迭代次数不足。

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

 cout << cdf(0.01) << endl;

如果我增加工作区的大小，则贝塞尔函数将不会计算。

我想知道是否有人可以让我深入了解我的问题。使用 x = 0.01 调用上面相应的 Mathematica 函数 F 返回 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 可能有点棘手；我对贝塞尔函数不是很熟悉，所以我帮不了你太多。

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