6

我想整合一个概率密度函数，(-\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返回一个错误，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;

如果我增加工作空间的大小，那么贝塞尔函数将不会计算。

我想知道是否有任何人可以让我了解我的问题。x = 0.01通过返回调用上述相应的 Mathematica 函数0.904384F。

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

密度

4

3 回答 3

1

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

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

计算出 K 可能有点棘手。我对贝塞尔函数不是很熟悉，所以在那里我帮不了你太多。

一般来说，我发现对于不明显的数值计算，最好的方法是在进行数值评估之前尽可能多地进行分析数学。（有点像自动对焦相机——把它靠近你想要的地方，然后让相机做剩下的事情。）

于 2012-05-28T03:56:49.013 回答
1

我在http://linux.math.tifr.res.in/manuals/html/gsl-ref-html/gsl-ref_16.html找到了关于这个 glsl 的完整描述 ，你可能会找到有用的信息。

由于我不是 GSL 专家，因此我没有从数学的角度关注您的问题，而是要提醒您有关浮点编程的一些关键方面。

您无法使用 IEEE 754 标准准确表示数字。MathLab 确实通过使用无限数字表示逻辑来隐藏事实，以便为您提供无错误的结果，这就是它与本机代码相比速度慢的原因。

我强烈推荐使用 FPU 参与科学微积分的任何人使用此链接：http: //docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

假设您喜欢那篇文章，我在上面的 GSL 链接上注意到了这一点：“如果错误界限太严格，例程将无法收敛”。

如果上限和下限之间的差异小于 double 的最小可表示值，即 std::numeric_limits::epsilon();，则您的界限可能过于严格。

另外请记住，从第二个链接开始，对于任何 C/C++ 编译器实现，默认舍入模式是“截断”，这会引入导致错误结果的细微微积分错误。我确实遇到了一个简单的 Liang Barsky 剪线器的问题，一阶！所以想象一下这一行的混乱：

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));

作为一般规则，在 C/C++ 中添加额外的变量来保存中间结果是明智的，这样您就可以逐步调试，然后看到任何舍入错误，您不应该尝试在任何本机编程中输入这样的表达式语言。不能比编译器更好地优化变量。

最后，作为一般规则，除非您对微积分的动态行为有信心，否则您应该先乘以再除。

祝你好运。

于 2012-05-28T02:46:32.030 回答
1

我没有尝试过 C++ 代码，但是通过检查 Mathematica 中的函数，它似乎在 mu 附近达到了极高的峰值，峰值的扩散由参数 lambda、alpha、beta 决定。

我要做的是对 pdf 进行初步搜索：查看 x=mu 的右侧和左侧，直到找到低于给定容差的第一个值。使用这些作为 cdf 的边界，而不是负无穷大。

伪代码如下：

x_mu
step = 0.000001
adaptive_step(y_value) -> returns a small step size if close to 0, and larger if far.

while (pdf_current > tolerance):
  step = adaptive_step(pdf_current)
  xtest = xtest - step
  pdf_current = pdf(xtest)

left_bound = xtest

//repeat for left bound

考虑到这个函数的峰值似乎有多紧密，收紧边界可能会为您节省大量当前浪费在计算零上的计算机时间。此外，您可以使用有界积分例程，而不是 -\infty,b 。

只是一个想法...

PS：Mathematica 给我 F[0.01, 3, 265, -5, 0] = 0.884505

于 2012-05-28T00:24:35.527 回答