StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

POQuadrature routines for probability densities
text
Body
copied!<p>I want to integrate a probability density function from <code>(-\infty, a]</code> because the cdf is not available in closed form. But I'm not sure how to do this in C++. </p> <p>This task is pretty simple in Mathematica; All I need to do is define the function, </p> <pre><code>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)) ]; </code></pre> <p>and then call the <code>NIntegrate</code> Routine to numerically integrate it. </p> <pre><code>F[x_, lambda_, alpha_, beta_, mu_] := NIntegrate[f[t, lambda, alpha, beta, mu], {t, -\[Infinity], x}] </code></pre> <p>Now I want to achieve the same thing in C++. I using the routine <code>gsl_integration_qagil</code> from the gsl numerics library. It is designed to integrate functions on the semi infinite intervals <code>(-\infty, a]</code> which is just what I want. But unfortunately I can't get it to work.</p> <p>This is the density function in C++,</p> <pre><code>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)); } </code></pre> <p>Then I try and integrate to obtain the cdf by calling the gsl routine. </p> <pre><code>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; } </code></pre> <p>However <code>gsl_integration_qagil</code> returns an error, <code>number of iterations was insufficient</code>. </p> <pre><code> double mu = 0.0f; double lambda = 3.0f; double alpha = 265.0f; double beta = -5.0f; cout << cdf(0.01) << endl; </code></pre> <p>If I increase the size of the workspace then the bessel function will not evaluate. </p> <p>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 <code>x = 0.01</code> returns <code>0.904384</code>. </p> <p>Could it be that the density is concentrated around a very small interval (i.e. outside of <code>[-0.05, 0.05]</code> the density is almost <code>0</code>, a plot is given below). If so what can be done about this. Thanks for reading. </p> <p><img src="https://i.stack.imgur.com/ZIgzK.png" alt="density"></p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload