Квадратурные подпрограммы для плотности вероятности

Я хочу проинтегрировать функцию плотности вероятности из (-\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_qagilиз библиотеки gsl numerics.Она предназначена для интегрирования функций на полубесконечных интервалах (-\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));

}  

Затем я пытаюсь интегрировать, чтобы получить cdf с помощью ca заполнение подпрограммы gsl.

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;

Если я увеличу размер рабочей области, функция Бесселя не будет оцениваться.

Мне было интересно, есть ли кто-нибудь, кто мог бы дать мне какое-либо представление о моей проблеме. Вызов соответствующей функции Mathematica F выше с x = 0,01возвращает 0,904384.

Может ли быть так, что плотность сосредоточена вокруг очень маленького интервала (то есть за пределами [-0,05, 0.05]плотность почти 0, график приведен ниже). Если да, то что можно сделать по этому поводу. Спасибо за чтение.