Quadrature routines for probability densities

Question

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

This task is quite simple in Mathematica; All I have to do is define a function,

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

and then call NIntegrate Routine for its numerical integration.

 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 am using the gsl_integration_qagil routine from the gsl library. It is designed to integrate functions at half-infinite intervals (-\infty, a] that I want. But, unfortunately, I cannot get it to work.

This is a 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 to integrate to get cdf by calling the gsl procedure.

 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, the bessel function will not be evaluated.

I was wondering if there is anyone who could let me understand my problem. Calling the corresponding Mathematica F function above with x = 0.01 returns 0.904384 .

Maybe the density is concentrated around a very small interval (that is, outside [-0.05, 0.05] density is almost 0 , the graph is shown below). If so, what can be done about this. Thanks for reading.

Answer 1

Re: integrating into +/- infinity:

I would use Mathematica to find an empirical estimate for | x is & mu; | → K, where K represents the "width" around the average value, and K is the function of alpha, beta and lambda - for example, F is less and approximately equal to (x- and mu;) ^{- 2} or ae ^{-b (x- & mu; ) 2} or something else. These functions have well-known integrals to infinity, for which you can evaluate empirically. Then you can integrate numerically into K and use the limited approximation to get from K to infinity.

Calculating K can be a little complicated; I am not very familiar with the functions of Bessel, so I can not help you.

In general, I found that for numerical calculation, which is not obvious, the best way is to do as much analytical mathematics as you can before performing a numerical evaluation. (Like a camera with autofocus - get closer to where you want, and then let the camera do the rest.)

Answer 2

I have not tried C ++ code, but, while checking the function in Mathematica, it seems extremely pointed around mu, with the spread of the peak determined by the parameters lambda, alpha, beta.

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

Following is the pseudo code:

 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 

Given how tightly this feature is achieved, tightening boundaries will probably save you a lot of time on a computer that is currently missing from the calculation of zeros. In addition, you can use the limited integration procedure, rather than - \ infty, b.

Just a thought ...

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

Answer 3

I found a full description of this glsl http://linux.math.tifr.res.in/manuals/html/gsl-ref-html/gsl-ref_16.html , you can find useful information.

Since I am not a GSL expert, I have not focused on your problem from a mathematical point of view, but I must remind you of some key aspect of floating point programming.

You cannot represent numbers accurately using the IEEE 754 standard. MathLab really hides this fact using infinite number logic to give you blunders without errors, so the reason for this is slow compared to native code.

I highly recommend this link to anyone involved in scientific calculus using the FPU: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Assuming you liked this article, I noticed this on the GSL link above: "Subroutines do not converge if error boundaries are too strict."

Your estimates may be too strict if the difference between the top and bottom is less than the minimum representable double value, that is, stand :: numeric_limits :: epsilon () ;.

In addition, remember that from the 2nd link for any implementation of the C / C ++ compiler, the default rounding mode is “truncated”, this leads to thin calculus errors that lead to incorrect results. I had a problem with a simple linear clipper Liang Barsky, 1st order! So imagine the mess on this line:

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

As a rule, in C / C ++ it is reasonable to add an additional variable while holding intermediate results, so you can debug step by step, and then see any rounding error, you should not try to enter an expression like this in any native programming language. It is not possible to optimize variables better than the compiler.

Finally, as a rule, you should multiply everything and then share if you are not sure about the dynamic behavior of your calculus.

Good luck.