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>
64 x = exp( ( lgamma(
a ) +
log( P ) ) /
a );
65 }
else if ( P > 0.95 ) {
66 x = -log1p( -P ) + lgamma(
a );
70 x = ( xg < -0.5 * sqrt(
a ) ) ?
a : sqrt(
a ) * xg +
a;
79 REAL8 lambda, dP, phi;
83 while ( keepgoing == 1 ) {
89 if ( dP == 0.0 ||
n++ > 32 ) {
94 lambda = dP /
fmax( 2.0 * fabs( dP /
x ), phi );
98 REAL8 step1 = -( (
a - 1.0 ) /
x - 1.0 ) * lambda * lambda / 4.0;
101 if ( fabs( step1 ) < 0.5 * fabs( step0 ) ) {
105 if (
x + step > 0 ) {
111 if ( fabs( step0 ) > 1
e-6 *
x || fabs( step0 * phi ) > 1
e-6 * P ) {
137 x = -
log(
Q ) + lgamma(
a );
138 }
else if (
Q > 0.95 ) {
139 x = exp( ( lgamma(
a ) + log1p( -
Q ) ) /
a );
143 x = ( xg < -0.5 * sqrt(
a ) ) ?
a : sqrt(
a ) * xg +
a;
152 REAL8 lambda, dQ, phi;
156 while ( keepgoing == 1 ) {
162 if ( dQ == 0.0 ||
n++ > 32 ) {
167 lambda = -dQ /
fmax( 2 * fabs( dQ /
x ), phi );
169 REAL8 step0 = lambda;
170 REAL8 step1 = -( (
a - 1 ) /
x - 1 ) * lambda * lambda / 4.0;
173 if ( fabs( step1 ) < 0.5 * fabs( step0 ) ) {
177 if (
x + step > 0 ) {
183 if ( fabs( step0 ) > 1
e-6 *
x ) {
202 if ( fabs( dP ) <= 0.425 ) {
207 pp = ( P < 0.5 ) ? P : 1.0 - P;
233 if ( fabs( dQ ) <= 0.425 ) {
239 pp = (
Q < 0.5 ) ?
Q : 1.0 -
Q;
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 ) );
392 return fmin( prob, 1.0 );
403 if ( fromzero == 1 ) {
405 REAL8 pk = gsl_ran_poisson_pdf( 0, halfdelta ) * gsl_cdf_chisq_P(
x, dof );
408 if ( (
REAL8 )counter < halfdelta ) {
413 P = gsl_ran_poisson_pdf( counter, halfdelta );
414 C = gsl_cdf_chisq_P(
x, dof + 2.0 * counter );
417 if ( !( ok == 1 && (
REAL8 )counter < halfdelta && dp >=
err * pk ) ) {
424 return fmin( prob, 1.0 );
437 return (
REAL4 )prob;
442 INT8 counter = (
INT8 )floor( halfdelta );
443 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
446 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) *
log(
x * 0.5 ) -
x * 0.5 - lgamma( dof * 0.5 + counter ) );
463 if ( fromzero == 1 ) {
468 if ( (
REAL8 )counter < halfdelta ) {
473 P = gsl_ran_poisson_pdf( counter, halfdelta );
478 if ( !( ok == 1 && (
REAL8 )counter < halfdelta && dp >=
err * pk ) ) {
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 ) );
513 return fmin( prob, 1.0 );
519 return fmin( prob, 1.0 );
532 return (
REAL4 )prob;
537 INT8 counter = (
INT8 )floor( halfdelta );
538 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
541 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) *
log(
x * 0.5 ) -
x * 0.5 - lgamma( dof * 0.5 + counter ) );
572 INT8 counter = (
INT8 )floor( halfdelta );
573 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
576 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) *
log(
x * 0.5 ) -
x * 0.5 - lgamma( dof * 0.5 + counter ) );
582 return fmin( prob, 1.0 );
588 return fmin( prob, 1.0 );
601 return (
REAL4 )prob;
606 INT8 counter = (
INT8 )floor( halfdelta );
607 REAL8 P = gsl_ran_poisson_pdf( counter, halfdelta );
610 REAL8 E = exp( ( dof * 0.5 + counter - 1.0 ) *
log(
x * 0.5 ) -
x * 0.5 - lgamma( dof * 0.5 + counter ) );
631 REAL8 dofint = 0.5 * dof - 1.0;
634 REAL8 logreal8min = -708.3964185322641;
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 );
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 ) );
643 if (
ul < logreal8min ) {
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;
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;
664 INT8 K = GSL_MAX_INT( 0, (
INT8 )floor( 0.5 * ( sqrt( dofint * dofint + 4.0 * dx ) - dofint ) ) );
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 );
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 );
678 while ( keep == 1 ) {
679 term *= ( dofint +
k ) *
k / dx;
681 if (
k <= 0 || term <=
epsval( sumK ) || keep != 1 ) {
689 while ( keep == 1 ) {
690 term /= ( dofint +
k ) *
k / dx;
692 if ( term <=
epsval( sumK ) || keep != 1 ) {
697 return 0.5 * exp( lntK +
log( sumK ) );
714 if (
delta == 0.0 ) {
715 return gsl_cdf_chisq_Pinv(
p, dof );
719 INT4 count_limit = 100;
731 while (
count < count_limit ) {
740 while ( worse == 0 ) {
741 if ( !( fabs( newF - pk ) > fabs( F - pk ) * ( 1.0 + crit ) && fabs( xk - xnew ) > crit * xk ) ) {
744 xnew = 0.5 * ( xnew + xk );
750 if ( !( fabs( h ) > crit * fabs( xk ) && fabs( h ) > crit ) ) {
757 fprintf( stderr,
"%s: Warning! ncx2inv(%g, %g, %g) failed to converge!\n",
__func__,
p, dof,
delta );
770 if (
delta == 0.0 ) {
771 return (
REAL4 )gsl_cdf_chisq_Pinv(
p, dof );
775 INT4 count_limit = 100;
787 while (
count < count_limit ) {
796 while ( worse == 0 ) {
797 if ( !( fabs( newF - pk ) > fabs( F - pk ) * ( 1.0 + crit ) && fabs( xk - xnew ) > crit * xk ) ) {
800 xnew = 0.5 * ( xnew + xk );
806 if ( !( fabs( h ) > crit * fabs( xk ) && fabs( h ) > crit ) ) {
813 fprintf( stderr,
"%s: Warning! ncx2inv_float() failed to converge!\n",
__func__ );
821 return mu - sigma * gsl_cdf_ugaussian_Qinv(
p );
828 REAL8 snr = ( value - dof ) / sqrt( 2.0 * dof );
#define __func__
log an I/O error, i.e.
static double double delta
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 binodeviance(REAL8 x, REAL8 np)
REAL8 matlab_gamma_inc(REAL8 x, REAL8 a, INT4 upper)
REAL8 twospect_tail(REAL8 r)
REAL8 unitGaussianSNR(REAL8 value, REAL8 dof)
REAL8 ncx2inv(REAL8 p, REAL8 dof, REAL8 delta)
Matlab's ncx2inv function.
INT4 cdf_gamma_Qinv(REAL8 *out, REAL8 Q, REAL8 a, REAL8 b)
REAL4 ncx2cdf_float(REAL4 x, REAL4 dof, REAL4 delta)
REAL8 cdf_gamma_P_usingmatlab(REAL8 x, REAL8 a, REAL8 b)
REAL4 ncx2inv_float(REAL8 p, REAL8 dof, REAL8 delta)
REAL8 ncx2cdf_withouttinyprob_withmatlabchi2cdf(REAL8 x, REAL8 dof, REAL8 delta)
REAL8 cdf_gamma_P(REAL8 x, REAL8 a, REAL8 b)
REAL8 cdf_chisq_Qinv(REAL8 Q, REAL8 nu)
REAL4 ncx2cdf_float_withouttinyprob_withmatlabchi2cdf(REAL4 x, REAL4 dof, REAL4 delta)
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...
REAL8 ncx2cdf_withouttinyprob(REAL8 x, REAL8 dof, REAL8 delta)
REAL8 cdf_gamma_Pinv(REAL8 P, REAL8 a, REAL8 b)
REAL8 cdf_ugaussian_Pinv(REAL8 P)
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...
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...
REAL8 ncx2pdf(REAL8 x, REAL8 dof, REAL8 delta)
REAL8 cdf_chisq_Pinv(REAL8 P, REAL8 nu)
INT4 cdf_ugaussian_Qinv(REAL8 *out, REAL8 Q)
REAL8 cdf_gamma_Q_usingmatlab(REAL8 x, REAL8 a, REAL8 b)
REAL8 norminv(REAL8 p, REAL8 mu, REAL8 sigma)
INT4 cdf_gamma_Q(REAL8 *out, REAL8 x, REAL8 a, REAL8 b)
REAL4 ncx2cdf_float_withouttinyprob(REAL4 x, REAL4 dof, REAL4 delta)
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_REAL4(assertion,...)
#define XLAL_CHECK_REAL8(assertion,...)
p
RUN ANALYSIS SCRIPT ###.