src/pgm_density.c
/*
 * Approximate the distribution function of PG(h, z).
 *
 * Note: The first term of the sum is evaluated before the loop to avoid
 * redundancy.
 */
Approximate the distribution function of PG(h, z).
Note: The first term of the sum is evaluated before the loop to avoid
redundancy.
Approximate the distribution function of PG ( h z ) . Note : The first term of the sum is evaluated before the loop to avoid redundancy .
double
pgm_polyagamma_cdf(double x, double h, double z)
{
    if (islessequal(x, 0.)) {
        return 0.;
    }
    else if (isinf(x)) {
        return 1.;
    }
    z = fabs(z);
    double c, zn;
    logcdf_func_t logcdf;
    struct cdf_args arg = {.a = h, .x = x, .z = z};
    if (z > 0.) {
        logcdf = invgauss_logcdf;
        c = h * log1p(exp(-z));
        arg.s2x = sqrt(x);
        zn = z;
    }
    else {
        logcdf = invgamma_logcdf;
        c = h * PGM_LOG2;
        arg.s2x = sqrt(2. * x);
        zn = 0.;
    }
    double sum = exp(c + logcdf(&arg));
    double sign = -1.;
    c -= pgm_lgamma(h);
    for (unsigned int n = 1; n < PGM_MAX_SERIES_TERMS; n++, sign = -sign, zn = z * n) {
        arg.a = 2 * n + h;
        double term = exp(c + pgm_lgamma(n + h) + logcdf(&arg) - pgm_lgamma(n + 1) - zn);
        double prev_sum = sum;
        sum += sign * term;
        if (PGM_ISCLOSE(sum, prev_sum, 0., DBL_EPSILON)) {
            break;
        }
    }
    return sum;
}
double pgm_polyagamma_cdf ( double x , double h , double z ) { if ( islessequal ( x , 0. ) ) { return 0. ; } else if ( isinf ( x ) ) { return 1. ; } z = fabs ( z ) ; double c , zn ; logcdf_func_t logcdf ; struct cdf_args arg = { . a = h , . x = x , . z = z } ; if ( z > 0. ) { logcdf = invgauss_logcdf ; c = h * log1p ( exp ( - z ) ) ; arg . s2x = sqrt ( x ) ; zn = z ; } else { logcdf = invgamma_logcdf ; c = h * PGM_LOG2 ; arg . s2x = sqrt ( 2. * x ) ; zn = 0. ; } double sum = exp ( c + logcdf ( & arg ) ) ; double sign = -1. ; c -= pgm_lgamma ( h ) ; for ( unsigned int n = 1 ; n < PGM_MAX_SERIES_TERMS ; n ++ , sign = - sign , zn = z * n ) { arg . a = 2 * n + h ; double term = exp ( c + pgm_lgamma ( n + h ) + logcdf ( & arg ) - pgm_lgamma ( n + 1 ) - zn ) ; double prev_sum = sum ; sum += sign * term ; if ( PGM_ISCLOSE ( sum , prev_sum , 0. , DBL_EPSILON ) ) { break ; } } return sum ; }
Approximate the distribution function of PG(h, z).
Approximate the distribution function of PG ( h z ) .