Rutinas de cuadratura para densidades de probabilidad

Question

Quiero integrar una función de densidad de probabilidad de (-\infty, a] porque el cdf no está disponible en forma cerrada. Pero no estoy seguro de cómo hacer esto en C++.

Esta tarea es bastante simple en Mathematica; Todo lo que necesito hacer es definir la función,

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

y luego llamar a la rutina NIntegrate para integrar numéricamente.

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

Ahora quiero lograr lo mismo en C++. Utilizo la rutina gsl_integration_qagil de la biblioteca de gsl numerics. Está diseñado para integrar funciones en los intervalos semi infinitos (-\infty, a] que es justo lo que quiero. Pero desafortunadamente no puedo hacer que funcione.

Ésta es la función de densidad en 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)); 

} 

entonces trato e integrar para obtener el CDF llamando a la rutina 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; 

} 

Sin embargo gsl_integration_qagil devuelve un 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; 

Si aumento el tamaño del espacio de trabajo, la función bessel no se evaluará.

Me preguntaba si había alguien que pudiera darme una idea de mi problema. Una llamada a la función F correspondiente de Mathematica con x = 0.01 devuelve 0.904384.

Podría ser que la densidad se concentre en un intervalo muy pequeño (es decir, fuera de [-0.05, 0.05], la densidad es casi 0, se proporciona un gráfico a continuación). Si es así, ¿qué se puede hacer al respecto? Gracias por leer.

Answer 1

Re: integrar a +/- infinito:

me gustaría utilizar Mathematica para encontrar una cota empírica para | x - μ | >> K, donde K representa el "ancho" alrededor de la media, y K es una función de alfa, beta y lambda; por ejemplo, F es menor que y aproximadamente igual aa (x- μ) ^-2 o ae ^{-b (x- μ)} o lo que sea. Estas funciones tienen integrales conocidas hasta el infinito, para lo cual se puede evaluar empíricamente. Luego puede integrar numéricamente a K y usar la aproximación limitada para pasar de K a infinito.

Averiguar K puede ser un poco complicado; No estoy muy familiarizado con las funciones de Bessel, así que no puedo ayudarte mucho allí.

En general, he encontrado que para el cálculo numérico que no es obvio, la mejor manera es hacer la mayor cantidad posible de matemáticas analíticas antes de hacer una evaluación numérica. (Más o menos como una cámara de enfoque automático: acércalo a donde quieras y deja que la cámara haga el resto.)

Answer 2

No he probado el código C++, pero por el control de la función en Mathematica, parece extremadamente alcanzó su punto máximo alrededor de mu, con la difusión del pico determinado por los parámetros lambda, alfa, beta .

Lo que haría sería hacer una búsqueda preliminar del pdf: mire a la derecha y a la izquierda de x = mu hasta que encuentre el primer valor debajo de una tolerancia dada. Úselos como los límites de su cdf, en lugar del infinito negativo.

pseudo código siguiente:

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 

Teniendo en cuenta lo bien alcanzó su punto máximo esta función parece ser, el endurecimiento de los límites probablemente ahorrará mucho tiempo de ordenador que está actualmente desperdicia el cálculo de ceros. Además, podría usar la rutina de integración limitada, en lugar de - \ infty, b.

Sólo un pensamiento ...

PS: Mathematica me da F [0,01, 3, 265, -5, 0] = 0.884505

Answer 3

He encontrado una descripción completa sobre este glsl allí http://linux.math.tifr.res.in/manuals/html/gsl-ref-html/gsl-ref_16.html, puede encontrar información útil.

Como no soy experto en GSL, no me enfoqué en su problema desde el punto de vista matemático, sino que debo recordarle algunos aspectos clave sobre la programación de punto flotante.

No puede representar números con precisión utilizando el estándar IEEE 754. MathLab esconde el hecho al usar una lógica de representación numérica infinita, para darle resultados sin errores en el enrutamiento, esta es la razón por la cual es lento en comparación con el código nativo.

Recomiendo encarecidamente este enlace para cualquier persona involucrada en el cálculo científico con un FPU: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Suponiendo que hayas disfrutado de este artículo, me he dado cuenta de esto en el enlace GSL arriba: "Las rutinas no podrán convergen si los límites de error son demasiado estrictos ".

Sus límites pueden ser demasiado estrictos si la diferencia entre la parte superior e inferior es menor que el valor mínimo representable de double, es decir std :: numeric_limits :: epsilon() ;.

Además, recuerde que, desde el segundo enlace, para cualquier implementación del compilador C/C++, el modo de redondeo predeterminado es "truncado", esto introduce errores de cálculo sutiles que dan lugar a resultados incorrectos. ¡Tuve el problema con una simple cortadora de línea Liang Barsky, primer orden! Así que imaginar el lío en esta línea:

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

Como regla general, es aconsejable en C/C++, para agregar variable adicional que sostiene los resultados intermedios, por lo que se puede depurar paso a paso, a continuación, ver ningún error de redondeo, no deberías intentar ingresar expresiones como esta en cualquier idioma nativo de programación. Uno no puede optimizar las variables mejor que un compilador.

Finalmente, como regla general, debe multiplicar todo y luego dividir, a menos que tenga confianza en el comportamiento dinámico de su cálculo.

Buena suerte.