src/pgm_saddle.c
/*
 * Sample from PG(h, z) using the Saddle approximation method.
 *
 *  Parameters
 *  ----------
 *  h : double
 *      The shape parameter of the distribution. The value must be a positive.
 *  z : double
 *      The exponential tilting parameter of the distributon.
 *  n : size_t
 *      The number of samples to generate.
 *  out: array of type double
 *      The array to place the generated samples. Only the first `n` elements
 *      will be populated.
 */
Sample from PG(h, z) using the Saddle approximation method.
Parameters

h : double
The shape parameter of the distribution. The value must be a positive.
z : double
The exponential tilting parameter of the distributon.
n : size_t
The number of samples to generate.
out: array of type double
The array to place the generated samples. Only the first `n` elements
will be populated.
Sample from PG ( h z ) using the Saddle approximation method . Parameters h : double The shape parameter of the distribution . The value must be a positive . z : double The exponential tilting parameter of the distributon . n : size_t The number of samples to generate . out : array of type double The array to place the generated samples . Only the first ` n ` elements will be populated .
void
random_polyagamma_saddle(bitgen_t* bitgen_state, double h, double z,
                         size_t n, double* out)
{
    double sqrt_rho, sqrt_rho_inv, hrho;
    float proposal_probability, p, q;
    parameter_t pr;
    set_sampling_parameters(&pr, h, 0.5 * fabs(z));
    sqrt_rho = sqrt(-2. * pr.left_tangent_slope);
    sqrt_rho_inv = 1. / sqrt_rho;
    p = expf(h * (0.5 / pr.xc + pr.left_tangent_intercept - sqrt_rho) +
             invgauss_logcdf(pr.xc, sqrt_rho_inv, h)) * pr.sqrt_alpha;
    hrho = -h * pr.right_tangent_slope;
    q = upper_incomplete_gamma(h, hrho * pr.xc, false) * pr.right_kernel_coef *
        expf(h * (pr.right_tangent_intercept - logf(hrho)));
    proposal_probability = p / (p + q);
    double mu2 = sqrt_rho_inv * sqrt_rho_inv;
    while (n--) {
        do {
            if (next_float(bitgen_state) < proposal_probability) {
                do {
                    double y = random_standard_normal(bitgen_state);
                    double w = sqrt_rho_inv + 0.5 * mu2 * y * y / h;
                    pr.x = w - sqrt(fabs(w * w - mu2));
                    if (next_double(bitgen_state) * (1. + pr.x * sqrt_rho) > 1.) {
                        pr.x = mu2 / pr.x;
                    }
                } while (pr.x >= pr.xc);
            }
            else {
                pr.x = random_left_bounded_gamma(bitgen_state, h, hrho, pr.xc);
            }
        } while (next_float(bitgen_state) * bounding_kernel(&pr) > saddle_point(&pr));
        out[n] = 0.25 * h * pr.x;
    }
}
void random_polyagamma_saddle ( bitgen_t * bitgen_state , double h , double z , size_t n , double * out ) { double sqrt_rho , sqrt_rho_inv , hrho ; float proposal_probability , p , q ; parameter_t pr ; set_sampling_parameters ( & pr , h , 0.5 * fabs ( z ) ) ; sqrt_rho = sqrt ( -2. * pr . left_tangent_slope ) ; sqrt_rho_inv = 1. / sqrt_rho ; p = expf ( h * ( 0.5 / pr . xc + pr . left_tangent_intercept - sqrt_rho ) + invgauss_logcdf ( pr . xc , sqrt_rho_inv , h ) ) * pr . sqrt_alpha ; hrho = - h * pr . right_tangent_slope ; q = upper_incomplete_gamma ( h , hrho * pr . xc , false ) * pr . right_kernel_coef * expf ( h * ( pr . right_tangent_intercept - logf ( hrho ) ) ) ; proposal_probability = p / ( p + q ) ; double mu2 = sqrt_rho_inv * sqrt_rho_inv ; while ( n -- ) { do { if ( next_float ( bitgen_state ) < proposal_probability ) { do { double y = random_standard_normal ( bitgen_state ) ; double w = sqrt_rho_inv + 0.5 * mu2 * y * y / h ; pr . x = w - sqrt ( fabs ( w * w - mu2 ) ) ; if ( next_double ( bitgen_state ) * ( 1. + pr . x * sqrt_rho ) > 1. ) { pr . x = mu2 / pr . x ; } } while ( pr . x >= pr . xc ) ; } else { pr . x = random_left_bounded_gamma ( bitgen_state , h , hrho , pr . xc ) ; } } while ( next_float ( bitgen_state ) * bounding_kernel ( & pr ) > saddle_point ( & pr ) ) ; out [ n ] = 0.25 * h * pr . x ; } }
Sample from PG(h, z) using the Saddle approximation method.
Sample from PG ( h z ) using the Saddle approximation method .