Routine di quadratura per densità di probabilità

Question

Desidero integrare una funzione di densità di probabilità da (-\infty, a] perché il cdf non è disponibile in forma chiusa. Ma non sono sicuro di come farlo in C++.

Questo compito è piuttosto semplice in Mathematica; Tutto quello che devi fare è definire la funzione,

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

e quindi chiamare il NIntegrate Routine di integrare numericamente esso.

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

Ora voglio ottenere la stessa cosa in C++. Io uso la routine gsl_integration_qagil dalla libreria numerica di gsl. È progettato per integrare le funzioni sugli intervalli semi infiniti (-\infty, a] che è proprio quello che voglio. Ma sfortunatamente non riesco a farlo funzionare.

Questa è la funzione di densità 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)); 

} 

Poi cerco di integrare per ottenere il CDF chiamando la routine di 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; 

} 

Tuttavia gsl_integration_qagil restituisce un errore, 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; 

Se si aumenta la dimensione dell'area di lavoro, la funzione di bessel non verrà valutata.

Mi chiedevo se c'era qualcuno che potesse darmi qualche idea del mio problema. Una chiamata alla corrispondente funzione Mathematica F sopra con x = 0.01 restituisce 0.904384.

Potrebbe essere che la densità si concentra su un intervallo molto piccolo (ad esempio all'esterno di [-0.05, 0.05] la densità è quasi 0, un grafico è riportato di seguito). In tal caso, cosa si può fare al riguardo. Grazie per aver letto.

Answer 1

Re: integrazione a +/- all'infinito:

avrei usato Mathematica per trovare un limite empirico per | x - μ | >> K, dove K rappresenta la "larghezza" attorno alla media, e K è una funzione di alfa, beta e lambda - per esempio F è minore di e approssimativamente uguale a (x- μ) ^-2 o ae ^{-b (x- μ)} o qualsiasi altra cosa. Queste funzioni hanno note integrali fino all'infinito, per le quali è possibile valutare empiricamente. Quindi puoi integrare numericamente a K e usare l'approssimazione limitata per ottenere da K all'infinito.

Capire K potrebbe essere un po 'complicato; Non ho molta familiarità con le funzioni di Bessel quindi non posso aiutarti molto lì.

In generale, ho trovato che per il calcolo numerico non è ovvio, il modo migliore è fare la matematica analitica più che puoi prima dello eseguire la valutazione numerica. (Un po 'come una fotocamera autofocus - avvicinala al punto desiderato, quindi lascia che la fotocamera faccia il resto.)

Answer 2

non hanno provato il codice C++, ma controllando la funzione in Mathematica, sembra estremamente picco intorno mu, con la diffusione del picco stabilito dai parametri lambda, alfa, beta .

Quello che vorrei fare sarebbe una ricerca preliminare del pdf: guarda a destra e a sinistra di x = mu finché non trovi il primo valore sotto una tolleranza data. Usa questi come limiti per il tuo cdf, invece dell'infinito negativo.

pseudo codice segue:

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 

Dato come ha raggiunto il picco strettamente questa funzione sembra essere, stringendo i limiti sarebbe probabilmente risparmiare un sacco di tempo al computer che è attualmente sprecato calcolo zeri. Inoltre, saresti in grado di utilizzare la routine di integrazione limitata, anziché - \ infty, b.

Solo un pensiero ...

PS: Mathematica mi dà F [0,01, 3, 265, -5, 0] = 0.884505

Answer 3

Ho trovato una descrizione completa su questo indirizzo http://linux.math.tifr.res.in/manuals/html/gsl-ref-html/gsl-ref_16.html, è possibile trovare informazioni utili.

Dato che non sono un esperto GSL, non mi sono concentrato sul problema dal punto di vista matematico, ma piuttosto devo ricordarvi alcuni aspetti chiave della programmazione in virgola mobile.

Non è possibile rappresentare in modo preciso i numeri utilizzando lo standard IEEE 754. MathLab nasconde il fatto utilizzando una logica di rappresentazione numerica infinita, al fine di darvi risultati senza errori, questo è il motivo per cui è lento rispetto al codice nativo.

Consiglierei vivamente questo link per chi si occupa di calcolo scientifico utilizzando una FPU: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Supponendo che avete goduto questo articolo, ho notato questo sul link GSL di cui sopra: "Le routine non riuscirà a convergere se i limiti di errore sono troppo stringenti ".

I limiti possono essere troppo stringenti se la differenza tra il superiore e il minore è inferiore al valore minimo rappresentabile del doppio, ovvero std :: numeric_limits :: epsilon() ;.

Inoltre, dal 2 ° collegamento, per qualsiasi implementazione del compilatore C/C++ la modalità di arrotondamento predefinita è "troncata", introducendo errori di calcolo sottili che sfiorano i risultati errati. Ho avuto il problema con un semplice tagliatore di Liang Barsky, 1 ° ordine! Quindi immaginate il caos in questa linea:

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

Come regola generale, è saggio in C/C++, aggiungere variabili aggiuntive tenendo risultati intermedi, in modo da poter eseguire il debug passo dopo passo, poi vedere alcun errore di arrotondamento, non dovresti provare a inserire espressioni come questa in qualsiasi lingua di programmazione nativa. Non è possibile ottimizzare le variabili meglio di un compilatore.

Infine, come regola generale, è necessario moltiplicare tutto e quindi dividere, a meno che non si abbia fiducia nel comportamento dinamico del calcolo.

Buona fortuna.