Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
cdfdist.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2014 Evan Goetz
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19
20//Based on GSL functions to determine chi-squared inversions
21//Some functions based from Matab 2012a functions, but optimized for TwoSpect analysis
22
23#include <math.h>
24#include <lal/LALConstants.h>
25#include <gsl/gsl_randist.h>
26#include <gsl/gsl_cdf.h>
27#include <gsl/gsl_sf_bessel.h>
28#include <gsl/gsl_math.h>
29#include "cdfdist.h"
30#include "TwoSpectSpecFunc.h"
31
33{
34 REAL8 val = cdf_gamma_Pinv( P, 0.5 * nu, 2.0 );
36 return val;
37}
39{
40 REAL8 val;
41 XLAL_CHECK_REAL8( cdf_gamma_Qinv( &val, Q, 0.5 * nu, 2.0 ) == XLAL_SUCCESS, XLAL_EFUNC );
42 return val;
43}
45{
46 REAL8 x;
47
48 XLAL_CHECK_REAL8( P < 1.0, XLAL_EFPOVRFLW, "Input P of 1.0 or larger returns infinity\n" );
49 if ( P == 0.0 ) {
50 return 0.0;
51 }
52
53 /* Consider, small, large and intermediate cases separately. The
54 boundaries at 0.05 and 0.95 have not been optimised, but seem ok
55 for an initial approximation.
56
57 BJG: These approximations aren't really valid, the relevant
58 criterion is P*gamma(a+1) < 1. Need to rework these routines and
59 use a single bisection style solver for all the inverse
60 functions.
61 */
62
63 if ( P < 0.05 ) {
64 x = exp( ( lgamma( a ) + log( P ) ) / a );
65 } else if ( P > 0.95 ) {
66 x = -log1p( -P ) + lgamma( a );
67 } else {
68 REAL8 xg = cdf_ugaussian_Pinv( P );
70 x = ( xg < -0.5 * sqrt( a ) ) ? a : sqrt( a ) * xg + a;
71 }
72
73 /* Use Lagrange's interpolation for E(x)/phi(x0) to work backwards
74 to an improved value of x (Abramowitz & Stegun, 3.6.6)
75
76 where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
77 */
78
79 REAL8 lambda, dP, phi;
80 UINT4 n = 0;
81
82 INT4 keepgoing = 1;
83 while ( keepgoing == 1 ) {
84 REAL8 val = cdf_gamma_P( x, a, 1.0 );
86 dP = P - val;
87 phi = ran_gamma_pdf( x, a, 1.0 );
88
89 if ( dP == 0.0 || n++ > 32 ) {
90 XLAL_CHECK_REAL8( fabs( dP ) <= sqrt( LAL_REAL4_EPS ) * P, XLAL_EFPINEXCT, "Inverse failed to converge\n" );
91 return b * x;
92 }
93
94 lambda = dP / fmax( 2.0 * fabs( dP / x ), phi );
95
96
97 REAL8 step0 = lambda;
98 REAL8 step1 = -( ( a - 1.0 ) / x - 1.0 ) * lambda * lambda / 4.0;
99
100 REAL8 step = step0;
101 if ( fabs( step1 ) < 0.5 * fabs( step0 ) ) {
102 step += step1;
103 }
104
105 if ( x + step > 0 ) {
106 x += step;
107 } else {
108 x *= 0.5;
109 }
110
111 if ( fabs( step0 ) > 1e-6 * x || fabs( step0 * phi ) > 1e-6 * P ) {
112 keepgoing = 1;
113 } else {
114 keepgoing = 0;
115 }
116 }
117
118 XLAL_CHECK_REAL8( fabs( dP ) <= sqrt( LAL_REAL4_EPS ) * P, XLAL_EFPINEXCT, "Inverse failed to converge\n" );
119 return b * x;
120
121}
123{
124 REAL8 x;
125
126 if ( Q == 1.0 ) {
127 *out = 0.0;
128 return XLAL_SUCCESS;
129 }
130 XLAL_CHECK( Q > 0.0 && Q < 1.0, XLAL_EFPOVRFLW, "Input P of 0.0 returns infinity\n" );
131
132 /* Consider, small, large and intermediate cases separately. The
133 boundaries at 0.05 and 0.95 have not been optimised, but seem ok
134 for an initial approximation. */
135
136 if ( Q < 0.05 ) {
137 x = -log( Q ) + lgamma( a );
138 } else if ( Q > 0.95 ) {
139 x = exp( ( lgamma( a ) + log1p( -Q ) ) / a );
140 } else {
141 REAL8 xg;
143 x = ( xg < -0.5 * sqrt( a ) ) ? a : sqrt( a ) * xg + a;
144 }
145
146 /* Use Lagrange's interpolation for E(x)/phi(x0) to work backwards
147 to an improved value of x (Abramowitz & Stegun, 3.6.6)
148
149 where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
150 */
151
152 REAL8 lambda, dQ, phi;
153 UINT4 n = 0;
154
155 BOOLEAN keepgoing = 1;
156 while ( keepgoing == 1 ) {
157 REAL8 val;
158 XLAL_CHECK( cdf_gamma_Q( &val, x, a, 1.0 ) == XLAL_SUCCESS, XLAL_EFUNC );
159 dQ = Q - val;
160 phi = ran_gamma_pdf( x, a, 1.0 );
161
162 if ( dQ == 0.0 || n++ > 32 ) {
163 *out = b * x;
164 return XLAL_SUCCESS;
165 }
166
167 lambda = -dQ / fmax( 2 * fabs( dQ / x ), phi );
168
169 REAL8 step0 = lambda;
170 REAL8 step1 = -( ( a - 1 ) / x - 1 ) * lambda * lambda / 4.0;
171
172 REAL8 step = step0;
173 if ( fabs( step1 ) < 0.5 * fabs( step0 ) ) {
174 step += step1;
175 }
176
177 if ( x + step > 0 ) {
178 x += step;
179 } else {
180 x /= 2.0;
181 }
182
183 if ( fabs( step0 ) > 1e-6 * x ) {
184 keepgoing = 1;
185 } else {
186 keepgoing = 0;
187 }
188 }
189
190 *out = b * x;
191 return XLAL_SUCCESS;
192}
194{
195 REAL8 r, x, pp;
196
197 REAL8 dP = P - 0.5;
198
199 XLAL_CHECK_REAL8( P != 1.0, XLAL_EFPOVRFLW );
200 XLAL_CHECK_REAL8( P != 0.0, XLAL_EFPOVRFLW );
201
202 if ( fabs( dP ) <= 0.425 ) {
203 x = twospect_small( dP );
204 return x;
205 }
206
207 pp = ( P < 0.5 ) ? P : 1.0 - P;
208
209 r = sqrt( -log( pp ) );
210
211 if ( r <= 5.0 ) {
213 } else {
214 x = twospect_tail( r );
215 }
216
217 if ( P < 0.5 ) {
218 return -x;
219 } else {
220 return x;
221 }
222
223}
225{
226 REAL8 r, x, pp;
227
228 REAL8 dQ = Q - 0.5;
229
230 XLAL_CHECK( Q != 1.0, XLAL_EFPOVRFLW );
231 XLAL_CHECK( Q != 0.0, XLAL_EFPOVRFLW );
232
233 if ( fabs( dQ ) <= 0.425 ) {
234 x = twospect_small( dQ );
235 *out = -x;
236 return XLAL_SUCCESS;
237 }
238
239 pp = ( Q < 0.5 ) ? Q : 1.0 - Q;
240
241 r = sqrt( -log( pp ) );
242
243 if ( r <= 5.0 ) {
245 } else {
246 x = twospect_tail( r );
247 }
248
249 if ( Q < 0.5 ) {
250 *out = x;
251 } else {
252 *out = -x;
253 }
254 return XLAL_SUCCESS;
255}
257{
258 REAL8 P;
259 REAL8 y = x / b;
260
261 if ( x <= 0.0 ) {
262 return 0.0;
263 }
264
265 if ( y > a ) {
266 REAL8 val;
268 P = 1.0 - val;
269 } else {
271 }
272
273 return P;
274}
276{
277
278 REAL8 P;
279 REAL8 y = x / b;
280
281 if ( x <= 0.0 ) {
282 return 0.0;
283 }
284
285 if ( y > a ) {
286 REAL8 val = matlab_gamma_inc( y, a, 1 );
287 P = 1.0 - val;
288 } else {
289 P = matlab_gamma_inc( y, a, 0 );
290 }
291
292 return P;
293}
295{
296 REAL8 Q;
297 REAL8 y = x / b;
298
299 if ( x <= 0.0 ) {
300 return 1.0;
301 }
302
303 if ( y < a ) {
304 REAL8 val;
306 Q = 1.0 - val;
307 } else {
309 }
310
311 *out = Q;
312 return XLAL_SUCCESS;
313}
315{
316
317 REAL8 Q;
318 REAL8 y = x / b;
319
320 if ( x <= 0.0 ) {
321 return 1.0;
322 }
323
324 if ( y < a ) {
325 REAL8 val = matlab_gamma_inc( y, a, 0 );
326 Q = 1.0 - val;
327 } else {
328 Q = matlab_gamma_inc( y, a, 1 );
329 }
330
331 return Q;
332}
333
334/**
335 * Compute the CDF P value at value x of a chi squared distribution with nu degrees of freedom
336 * Rougly REAL4 precision
337 * \param [in] x CDF value at value x
338 * \param [in] nu Number of degrees of freedom
339 * \return CDF value
340 */
342{
343 REAL8 val = cdf_gamma_P( x, 0.5 * nu, 2.0 );
345 return val;
346} /* twospect_cdf_chisq_P() */
347
348
349/**
350 * Compute the CDF P value at value x of a chi squared distrubution with nu degrees of freedom using the Matlab-based function
351 * \param [in] x CDF value at value x
352 * \param [in] nu Number of degrees of freedom
353 * \return CDF value
354 */
356{
357 REAL8 val = cdf_gamma_P_usingmatlab( x, 0.5 * nu, 2.0 );
359 return val;
360} /* matlab_cdf_chisq_P() */
361
362
363/**
364 * Matlab's version of the non-central chi-squared CDF with nu degrees of freedom and non-centrality delta at value x
365 * \param [in] x Value at which to compute the CDF
366 * \param [in] dof Number of degrees of freedom
367 * \param [in] delta Non-centrality parameter
368 * \return CDF value
369 */
371{
372
373 REAL8 prob = 0.0;
374
375 //Fail for bad inputs or return 0 if x<=0
376 XLAL_CHECK_REAL8( dof >= 0.0 && delta >= 0.0, XLAL_EINVAL );
377 if ( x <= 0.0 ) {
378 return prob;
379 }
380
382 REAL8 halfdelta = 0.5 * delta;
383 INT8 counter = ( INT8 )floor( halfdelta );
384 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
385 REAL8 C = gsl_cdf_chisq_P( x, dof + 2.0 * counter );
386 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) * log( x * 0.5 ) - x * 0.5 - lgamma( dof * 0.5 + counter ) );
387
388 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 0 );
390 counter--;
391 if ( counter < 0 ) {
392 return fmin( prob, 1.0 );
393 }
394
395 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 1 );
397
398 //This part computes small probabilities
399 INT4 fromzero = 0;
400 if ( prob == 0.0 ) {
401 fromzero = 1;
402 }
403 if ( fromzero == 1 ) {
404 counter = 0;
405 REAL8 pk = gsl_ran_poisson_pdf( 0, halfdelta ) * gsl_cdf_chisq_P( x, dof );
406 REAL8 dp = 0.0;
407 INT4 ok = 0;
408 if ( ( REAL8 )counter < halfdelta ) {
409 ok = 1;
410 }
411 while ( ok == 1 ) {
412 counter++;
413 P = gsl_ran_poisson_pdf( counter, halfdelta );
414 C = gsl_cdf_chisq_P( x, dof + 2.0 * counter );
415 dp = P * C;
416 pk += dp;
417 if ( !( ok == 1 && ( REAL8 )counter < halfdelta && dp >= err * pk ) ) {
418 ok = 0;
419 }
420 }
421 prob = pk;
422 }
423
424 return fmin( prob, 1.0 );
425
426} /* ncx2cdf() */
427
428//Matlab's non-central chi square CDF up to REAL4 precision
430{
431
432 REAL8 prob = 0.0;
433
434 //Fail for bad inputs or return 0 if x<=0
435 XLAL_CHECK_REAL4( dof >= 0.0 && delta >= 0.0, XLAL_EINVAL );
436 if ( x <= 0.0 ) {
437 return ( REAL4 )prob;
438 }
439
441 REAL8 halfdelta = 0.5 * delta;
442 INT8 counter = ( INT8 )floor( halfdelta );
443 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
444 REAL8 C = twospect_cdf_chisq_P( ( REAL8 )x, ( REAL8 )( dof + 2.0 * counter ) );
446 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) * log( x * 0.5 ) - x * 0.5 - lgamma( dof * 0.5 + counter ) );
447
448 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 0 );
450 counter--;
451 if ( counter < 0 ) {
452 return ( REAL4 )fmin( prob, 1.0 );
453 }
454
455 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 1 );
457
458 //This part computes small probabilities
459 INT4 fromzero = 0;
460 if ( prob == 0.0 ) {
461 fromzero = 1;
462 }
463 if ( fromzero == 1 ) {
464 counter = 0;
465 REAL8 pk = gsl_ran_poisson_pdf( 0, halfdelta ) * twospect_cdf_chisq_P( x, dof );
466 REAL8 dp = 0.0;
467 INT4 ok = 0;
468 if ( ( REAL8 )counter < halfdelta ) {
469 ok = 1;
470 }
471 while ( ok == 1 ) {
472 counter++;
473 P = gsl_ran_poisson_pdf( counter, halfdelta );
474 C = twospect_cdf_chisq_P( x, dof + 2.0 * counter );
476 dp = P * C;
477 pk += dp;
478 if ( !( ok == 1 && ( REAL8 )counter < halfdelta && dp >= err * pk ) ) {
479 ok = 0;
480 }
481 }
482 prob = pk;
483 }
484
485 return ( REAL4 )fmin( prob, 1.0 );
486
487} /* ncx2cdf_float() */
488
489//Matlab's non-central chi-square tries to compute very small probabilities. We don't normally need this,
490//so this function leaves out the last part to compute small probabilities.
492{
493
494 REAL8 prob = 0.0;
495
496 //Fail for bad inputs or return 0 if x<=0
497 XLAL_CHECK_REAL8( dof >= 0.0 && delta >= 0.0, XLAL_EINVAL );
498 if ( x <= 0.0 ) {
499 return prob;
500 }
501
503 REAL8 halfdelta = 0.5 * delta;
504 INT8 counter = ( INT8 )floor( halfdelta );
505 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
506 REAL8 C = gsl_cdf_chisq_P( x, dof + 2.0 * counter );
507 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) * log( x * 0.5 ) - x * 0.5 - lgamma( dof * 0.5 + counter ) );
508
509 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 0 );
511 counter--;
512 if ( counter < 0 ) {
513 return fmin( prob, 1.0 );
514 }
515
516 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 1 );
518
519 return fmin( prob, 1.0 );
520
521} /* ncx2cdf_withouttinyprob() */
522
523//Without small probabilities up to REAL4 precision
525{
526
527 REAL8 prob = 0.0;
528
529 //Fail for bad inputs or return 0 if x<=0
530 XLAL_CHECK_REAL4( dof >= 0.0 && delta >= 0.0, XLAL_EINVAL );
531 if ( x <= 0.0 ) {
532 return ( REAL4 )prob;
533 }
534
536 REAL8 halfdelta = 0.5 * delta;
537 INT8 counter = ( INT8 )floor( halfdelta );
538 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
539 REAL8 C = twospect_cdf_chisq_P( ( REAL8 )x, ( REAL8 )( dof + 2.0 * counter ) );
541 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) * log( x * 0.5 ) - x * 0.5 - lgamma( dof * 0.5 + counter ) );
542
543 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 0 );
545 counter--;
546 if ( counter < 0 ) {
547 return ( REAL4 )fmin( prob, 1.0 );
548 }
549
550 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 1 );
552
553 return ( REAL4 )fmin( prob, 1.0 );
554
555} /* ncx2cdf_float_withouttinyprob() */
556
557
558//This is ncx2cdf function like in Matlab, but using the Matlab version of the central chi square calculation instead of the GSL version
560{
561
562 REAL8 prob = 0.0;
563
564 //Fail for bad inputs or return 0 if x<=0
565 XLAL_CHECK_REAL8( dof >= 0.0 && delta >= 0.0, XLAL_EINVAL );
566 if ( x <= 0.0 ) {
567 return prob;
568 }
569
571 REAL8 halfdelta = 0.5 * delta;
572 INT8 counter = ( INT8 )floor( halfdelta );
573 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
574 REAL8 C = matlab_cdf_chisq_P( x, dof + 2.0 * counter ); //Matlab chi square cdf calculation
576 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) * log( x * 0.5 ) - x * 0.5 - lgamma( dof * 0.5 + counter ) );
577
578 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 0 );
580 counter--;
581 if ( counter < 0 ) {
582 return fmin( prob, 1.0 );
583 }
584
585 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 1 );
587
588 return fmin( prob, 1.0 );
589
590} /* ncx2cdf_withouttinyprob_withmatlabchi2cdf() */
591
592//This is ncx2cdf function like in Matlab, but using the Matlab version of the central chi square calculation instead of the GSL version; up to REAL4 precision
594{
595
596 REAL8 prob = 0.0;
597
598 //Fail for bad inputs or return 0 if x<=0
599 XLAL_CHECK_REAL4( dof >= 0.0 && delta >= 0.0, XLAL_EINVAL );
600 if ( x <= 0.0 ) {
601 return ( REAL4 )prob;
602 }
603
605 REAL8 halfdelta = 0.5 * delta;
606 INT8 counter = ( INT8 )floor( halfdelta );
607 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
608 REAL8 C = matlab_cdf_chisq_P( ( REAL8 )x, ( REAL8 )( dof + 2.0 * counter ) ); //Matlab chi2cdf
610 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) * log( x * 0.5 ) - x * 0.5 - lgamma( dof * 0.5 + counter ) );
611
612 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 0 );
614 counter--;
615 if ( counter < 0 ) {
616 return ( REAL4 )fmin( prob, 1.0 );
617 }
618
619 sumseries_eg( &prob, P, C, E, counter, x, dof, halfdelta, err, 1 );
621
622 return ( REAL4 )fmin( prob, 1.0 );
623
624} /* ncx2cdf_float_withouttinyprob_withmatlabchi2cdf() */
625
626
627//Like Matlabs ncx2pdf
629{
630
631 REAL8 dofint = 0.5 * dof - 1.0;
632 REAL8 x1 = sqrt( x );
633 REAL8 delta1 = sqrt( delta );
634 REAL8 logreal8min = -708.3964185322641;
635
636 REAL8 ul = 0.0;
637 if ( dofint <= -0.5 ) {
638 ul = -0.5 * ( delta + x ) + 0.5 * x1 * delta1 / ( dofint + 1.0 ) + dofint * ( log( x ) - LAL_LN2 ) - LAL_LN2 - lgamma( dofint + 1.0 );
639 } else {
640 ul = -0.5 * ( delta1 - x1 ) * ( delta1 - x1 ) + dofint * ( log( x ) - LAL_LN2 ) - LAL_LN2 - lgamma( dofint + 1.0 ) + ( dofint + 0.5 ) * log( ( dofint + 0.5 ) / ( x1 * delta1 + dofint + 0.5 ) );
641 }
642
643 if ( ul < logreal8min ) {
644 return 0.0;
645 }
646
647 //Scaled Bessel function?
648 gsl_sf_result sbes = {0, 0};
649 INT4 status = gsl_sf_bessel_Inu_scaled_e( dofint, delta1 * x1, &sbes );
650 if ( status == GSL_SUCCESS && sbes.val > 0.0 ) {
651 return exp( -LAL_LN2 - 0.5 * ( x1 - delta1 ) * ( x1 - delta1 ) + dofint * log( x1 / delta1 ) ) * sbes.val;
652 }
653
654 //Bessel function without scaling?
655 gsl_sf_result bes;
656 status = gsl_sf_bessel_Inu_e( dofint, delta1 * x1, &bes );
657 if ( status == GSL_SUCCESS && bes.val > 0.0 ) {
658 return exp( -LAL_LN2 - 0.5 * ( x + delta ) + dofint * log( x1 / delta1 ) ) * bes.val;
659 }
660
661 //Okay, now recursion
662 REAL8 lnsr2pi = log( sqrt( LAL_TWOPI ) );
663 REAL8 dx = delta * x * 0.25;
664 INT8 K = GSL_MAX_INT( 0, ( INT8 )floor( 0.5 * ( sqrt( dofint * dofint + 4.0 * dx ) - dofint ) ) );
665 REAL8 lntK = 0.0;
666 if ( K == 0 ) {
667 lntK = -lnsr2pi - 0.5 * ( delta + log( dofint ) ) - ( lgamma( dofint + 1 ) - 0.5 * log( LAL_TWOPI * dofint ) + dofint * log( dofint ) - dofint ) - binodeviance( dofint, 0.5 * x );
668 } else {
669 lntK = -2.0 * lnsr2pi - 0.5 * ( log( K ) + log( dofint + K ) ) - ( lgamma( K + 1 ) - 0.5 * log( LAL_TWOPI * K ) + K * log( K ) - K ) - ( lgamma( dofint + K + 1 ) - 0.5 * log( LAL_TWOPI * ( dofint + K ) ) + ( dofint + K ) * log( dofint + K ) - ( dofint + K ) ) - binodeviance( K, 0.5 * delta ) - binodeviance( dofint + K, 0.5 * x );
670 }
671 REAL8 sumK = 1.0;
672 INT4 keep = 0;
673 if ( K > 0 ) {
674 keep = 1;
675 }
676 REAL8 term = 1.0;
677 REAL8 k = K;
678 while ( keep == 1 ) {
679 term *= ( dofint + k ) * k / dx;
680 sumK += term;
681 if ( k <= 0 || term <= epsval( sumK ) || keep != 1 ) {
682 keep = 0;
683 }
684 k--;
685 }
686 keep = 1;
687 term = 1.0;
688 k = K + 1;
689 while ( keep == 1 ) {
690 term /= ( dofint + k ) * k / dx;
691 sumK += term;
692 if ( term <= epsval( sumK ) || keep != 1 ) {
693 keep = 0;
694 }
695 k++;
696 }
697 return 0.5 * exp( lntK + log( sumK ) );
698
699} /* ncx2pdf() */
700
701/**
702 * Matlab's ncx2inv function
703 * \param [in] p CDF P value from which to compute the inversion
704 * \param [in] dof Number of degrees of freedom
705 * \param [in] delta Non-centrality parameter
706 * \return The x value that corresponds to the P value
707 */
709{
710
711 //Fail if bad input
713
714 if ( delta == 0.0 ) {
715 return gsl_cdf_chisq_Pinv( p, dof );
716 }
717
718 REAL8 pk = p;
719 INT4 count_limit = 100;
720 INT4 count = 0;
721 REAL8 crit = sqrt( LAL_REAL8_EPS );
722 REAL8 mn = dof + delta;
723 REAL8 variance = 2.0 * ( dof + 2.0 * delta );
724 REAL8 temp = log( variance + mn * mn );
725 REAL8 mu = 2.0 * log( mn ) - 0.5 * temp;
726 REAL8 sigma = -2.0 * log( mn ) + temp;
727 REAL8 xk = exp( norminv( pk, mu, sigma ) );
728 REAL8 h = 0.0;
729 REAL8 F = ncx2cdf( xk, dof, delta );
731 while ( count < count_limit ) {
732 count++;
733 REAL8 f = ncx2pdf( xk, dof, delta );
735 h = ( F - pk ) / f;
736 REAL8 xnew = fmax( 0.2 * xk, fmin( 5.0 * xk, xk - h ) );
737 REAL8 newF = ncx2cdf( xnew, dof, delta );
739 INT4 worse = 0;
740 while ( worse == 0 ) {
741 if ( !( fabs( newF - pk ) > fabs( F - pk ) * ( 1.0 + crit ) && fabs( xk - xnew ) > crit * xk ) ) {
742 worse = 1;
743 } else {
744 xnew = 0.5 * ( xnew + xk );
745 newF = ncx2cdf( xnew, dof, delta );
747 }
748 }
749 h = xk - xnew;
750 if ( !( fabs( h ) > crit * fabs( xk ) && fabs( h ) > crit ) ) {
751 return xk;
752 }
753 xk = xnew;
754 F = newF;
755 }
756
757 fprintf( stderr, "%s: Warning! ncx2inv(%g, %g, %g) failed to converge!\n", __func__, p, dof, delta );
758 return xk;
759
760} /* ncx2inv() */
761
762
763//Matlab's ncx2inv() function to REAL4 precision
765{
766
767 //Fail if bad input
769
770 if ( delta == 0.0 ) {
771 return ( REAL4 )gsl_cdf_chisq_Pinv( p, dof );
772 }
773
774 REAL8 pk = p;
775 INT4 count_limit = 100;
776 INT4 count = 0;
777 REAL8 crit = sqrt( LAL_REAL4_EPS );
778 REAL8 mn = dof + delta;
779 REAL8 variance = 2.0 * ( dof + 2.0 * delta );
780 REAL8 temp = log( variance + mn * mn );
781 REAL8 mu = 2.0 * log( mn ) - 0.5 * temp;
782 REAL8 sigma = -2.0 * log( mn ) + temp;
783 REAL8 xk = exp( norminv( pk, mu, sigma ) );
784 REAL8 h = 0.0;
787 while ( count < count_limit ) {
788 count++;
789 REAL8 f = ncx2pdf( xk, dof, delta );
791 h = ( F - pk ) / f;
792 REAL8 xnew = fmax( 0.2 * xk, fmin( 5.0 * xk, xk - h ) );
795 INT4 worse = 0;
796 while ( worse == 0 ) {
797 if ( !( fabs( newF - pk ) > fabs( F - pk ) * ( 1.0 + crit ) && fabs( xk - xnew ) > crit * xk ) ) {
798 worse = 1;
799 } else {
800 xnew = 0.5 * ( xnew + xk );
803 }
804 }
805 h = xk - xnew;
806 if ( !( fabs( h ) > crit * fabs( xk ) && fabs( h ) > crit ) ) {
807 return xk;
808 }
809 xk = xnew;
810 F = newF;
811 }
812
813 fprintf( stderr, "%s: Warning! ncx2inv_float() failed to converge!\n", __func__ );
814 return xk;
815
816} /* ncx2inv_float() */
817
818//Matlab's norminv function
820{
821 return mu - sigma * gsl_cdf_ugaussian_Qinv( p );
822} /* norminv() */
823
824
825//For the normal distribution, what is the SNR of a given value
827{
828 REAL8 snr = ( value - dof ) / sqrt( 2.0 * dof );
829 return snr;
830} /* unitGaussianSNR() */
#define __func__
log an I/O error, i.e.
int k
static double double delta
#define fprintf
void sumseries_eg(REAL8 *computedprob, REAL8 P, REAL8 C, REAL8 E, INT8 counter, REAL8 x, REAL8 dof, REAL8 halfdelta, REAL8 err, INT4 countdown)
REAL8 ran_gamma_pdf(REAL8 x, REAL8 a, REAL8 b)
REAL8 twospect_intermediate(REAL8 r)
REAL8 twospect_small(REAL8 q)
INT4 sf_gamma_inc_Q(REAL8 *out, REAL8 a, REAL8 x)
INT4 sf_gamma_inc_P(REAL8 *out, REAL8 a, REAL8 x)
REAL8 epsval(REAL8 val)
REAL8 binodeviance(REAL8 x, REAL8 np)
REAL8 matlab_gamma_inc(REAL8 x, REAL8 a, INT4 upper)
REAL8 twospect_tail(REAL8 r)
double e
REAL8 unitGaussianSNR(REAL8 value, REAL8 dof)
Definition: cdfdist.c:826
REAL8 ncx2inv(REAL8 p, REAL8 dof, REAL8 delta)
Matlab's ncx2inv function.
Definition: cdfdist.c:708
INT4 cdf_gamma_Qinv(REAL8 *out, REAL8 Q, REAL8 a, REAL8 b)
Definition: cdfdist.c:122
REAL4 ncx2cdf_float(REAL4 x, REAL4 dof, REAL4 delta)
Definition: cdfdist.c:429
REAL8 cdf_gamma_P_usingmatlab(REAL8 x, REAL8 a, REAL8 b)
Definition: cdfdist.c:275
REAL4 ncx2inv_float(REAL8 p, REAL8 dof, REAL8 delta)
Definition: cdfdist.c:764
REAL8 ncx2cdf_withouttinyprob_withmatlabchi2cdf(REAL8 x, REAL8 dof, REAL8 delta)
Definition: cdfdist.c:559
REAL8 cdf_gamma_P(REAL8 x, REAL8 a, REAL8 b)
Definition: cdfdist.c:256
REAL8 cdf_chisq_Qinv(REAL8 Q, REAL8 nu)
Definition: cdfdist.c:38
REAL4 ncx2cdf_float_withouttinyprob_withmatlabchi2cdf(REAL4 x, REAL4 dof, REAL4 delta)
Definition: cdfdist.c:593
REAL8 matlab_cdf_chisq_P(REAL8 x, REAL8 nu)
Compute the CDF P value at value x of a chi squared distrubution with nu degrees of freedom using the...
Definition: cdfdist.c:355
REAL8 ncx2cdf_withouttinyprob(REAL8 x, REAL8 dof, REAL8 delta)
Definition: cdfdist.c:491
REAL8 cdf_gamma_Pinv(REAL8 P, REAL8 a, REAL8 b)
Definition: cdfdist.c:44
REAL8 cdf_ugaussian_Pinv(REAL8 P)
Definition: cdfdist.c:193
REAL8 ncx2cdf(REAL8 x, REAL8 dof, REAL8 delta)
Matlab's version of the non-central chi-squared CDF with nu degrees of freedom and non-centrality del...
Definition: cdfdist.c:370
REAL8 twospect_cdf_chisq_P(REAL8 x, REAL8 nu)
Compute the CDF P value at value x of a chi squared distribution with nu degrees of freedom Rougly RE...
Definition: cdfdist.c:341
REAL8 ncx2pdf(REAL8 x, REAL8 dof, REAL8 delta)
Definition: cdfdist.c:628
REAL8 cdf_chisq_Pinv(REAL8 P, REAL8 nu)
Definition: cdfdist.c:32
INT4 cdf_ugaussian_Qinv(REAL8 *out, REAL8 Q)
Definition: cdfdist.c:224
REAL8 cdf_gamma_Q_usingmatlab(REAL8 x, REAL8 a, REAL8 b)
Definition: cdfdist.c:314
REAL8 norminv(REAL8 p, REAL8 mu, REAL8 sigma)
Definition: cdfdist.c:819
INT4 cdf_gamma_Q(REAL8 *out, REAL8 x, REAL8 a, REAL8 b)
Definition: cdfdist.c:294
REAL4 ncx2cdf_float_withouttinyprob(REAL4 x, REAL4 dof, REAL4 delta)
Definition: cdfdist.c:524
const double pp
const double Q
#define LAL_LN2
#define LAL_REAL8_EPS
#define LAL_TWOPI
#define LAL_REAL4_EPS
unsigned char BOOLEAN
double REAL8
int64_t INT8
uint32_t UINT4
int32_t INT4
float REAL4
static const INT4 r
static const INT4 a
#define xlalErrno
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_REAL4(assertion,...)
#define XLAL_CHECK_REAL8(assertion,...)
XLAL_EFPINEXCT
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EINVAL
XLAL_EFPOVRFLW
long long count
Definition: hwinject.c:371
list y
list mu
temp
n
out