Take the 2-minute tour ×

I want to integrate a probability density function from (-\infty, a] because the cdf is not available in closed form. But I'm not sure how to do this in C++.

This task is pretty simple in Mathematica; All I need to do is define the function,

f[x_, lambda_, alpha_, beta_, mu_] := 
     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))

and then call the NIntegrate Routine to numerically integrate it.

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

Now I want to achieve the same thing in C++. I using the routine gsl_integration_qagil from the gsl numerics library. It is designed to integrate functions on the semi infinite intervals (-\infty, a] which is just what I want. But unfortunately I can't get it to work.

This is the density function in 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));


Then I try and integrate to obtain the cdf by calling the gsl routine.

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;


However gsl_integration_qagil returns an error, 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;

If I increase the size of the workspace then the bessel function will not evaluate.

I was wondering if there was anyone that could give me any insight to my problem. A call to the corresponding Mathematica function F above with x = 0.01 returns 0.904384.

Could it be that the density is concentrated around a very small interval (i.e. outside of [-0.05, 0.05] the density is almost 0, a plot is given below). If so what can be done about this. Thanks for reading.


share|improve this question
This function looks symmetric, meaning that cdf(0) = 1/2. Remember that the cdf evaluated at x is the same as the integral from 0 to x plus the cdf evaluated at 0. Of course, I'm going from the shape of the graph, it may not actually be exactly symmetric. –  Ben Voigt May 28 '12 at 4:00
Hi @Ben, I wish that was the case! –  aukie May 28 '12 at 4:20

3 Answers 3

up vote 1 down vote accepted

Re: integrating to +/- infinity:

I would use Mathematica to find an empirical bound for |x - μ| >> K, where K represents the "width" around the mean, and K is a function of alpha, beta, and lambda -- for example F is less than and approximately equal to a(x-μ)-2 or ae-b(x-μ)2 or whatever. These functions have known integrals out to infinity, for which you can evaluate empirically. Then you can integrate numerically out to K, and use the bounded approximation to get from K to infinity.

Figuring out K may be a bit tricky; I'm not very familiar with Bessel functions so I can't help you much there.

In general, I've found that for numerical calculation that's not obvious, the best way is to do as much analytical math as you can before you do numerical evaluation. (Kind of like an autofocus camera -- get it close to where you want, then let the camera do the rest.)

share|improve this answer
Hi @Jason, given that I know the asymptotic behaviour of the itegral I may have to adopt this kind of approach. But I just don't get how mathematica worked out the integral. Maybe it is doing something similar behind the scenes. –  aukie May 28 '12 at 4:22
Do you have a copy of Numerical Recipes? I think maybe they discuss integrals that extend to infinity, not sure how they handle it. –  Jason S May 28 '12 at 12:25
Leafing through it now! –  aukie May 29 '12 at 19:41

I found a complete description on this glsl there http://linux.math.tifr.res.in/manuals/html/gsl-ref-html/gsl-ref_16.html, you may find usefull informations.

Since I'm not GSL expert I did not focus on your problem from the math point of view, but rather I've to remind you some key aspect about floating point programming.

You can't accurately represent numbers using IEEE 754 standard. MathLab do hide the fact by using an infinite number representation logic, in order to give you rouding error-free results , this is the reason why it's slow compared to native code.

I strongly recommand this link for anyone involved in scientific calculus using a FPU: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Assuming you've enjoyed that article, I've noticed this on the GSL link above: "The routines will fail to converge if the error bounds are too stringent".

Your bounds may be too stringent if the difference between the upper and the lower is less than the minimum representable value of double, that is std::numeric_limits::epsilon();.

In addition remember, from the 2nd link, for any C/C++ compiler implementation the default rounding mode is "truncate", this introduce subtle calculus errors leeding to the wrong results. I did have the problem with a simple Liang Barsky line clipper, 1st order ! So imagine the mess in this line:


As a general rule, it is wise in C/C++, to add additional variable holding intermediate results, so you can debug step by step, then see any rounding error, you shouldn't try to input expression like this one in any native programing langage. One can't optimize variables better than a compiler.

Finally, as a general rule, you should multiply everything then divide, unless you are confident about the dynamic behavior of your calculus.

Good luck.

share|improve this answer
Hi @Gold, I hate round off error! To paraphrase Moore, "rigour is lost when we need it most". He's referring to mathematicians implementing they're algorithms using floating point arithmetic. Hence the advent of interval arithmetic. I for one look forward to the day this becomes a hardware feature and we can say good bye to round off! (well not quite but at least have better understanding of it) I knew of that article before, but u've motivated me to go read it! Thanks for your input –  aukie May 28 '12 at 4:29

I haven't tried the C++ code, but by checking out the function in Mathematica, it does seem extremely peaked around mu, with the spread of the peak determined by the parameters lambda,alpha,beta.

What I would do would be to do a preliminary search of the pdf: look to the right and to the left of x=mu until you find the first value below a given tolerance. Use these as the bounds for your cdf, instead of negative infinity.

Pseudo code follows:

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

Given how tightly peaked this function seems to be, tightening the bounds would probably save you a lot of computer time that's currently wasted calculating zeros. Also, you'd be able to use the bounded integration routine, rather than -\infty,b .

Just a thought...

PS: Mathematica gives me F[0.01, 3, 265, -5, 0] = 0.884505

share|improve this answer
Hi @Adriano, Thanks for your input. I might just have to approach my problem in this kind of manner. –  aukie May 28 '12 at 4:17

Your Answer


