Function normal_pdf 

pub fn normal_pdf(x: f64) -> f64

Standard normal PDF phi(x).

The squared argument is carried exactly (see square_residual); without that, exp(-½·fl(x*x)) degrades like x²·ε/2 and reaches 5.7e-14 relative before φ underflows, against the 3.3e-16 it holds with.

Examples found in repository

examples/special_audit_dump.rs (line 78)

fn main() {
    // ---- Bessel channels -------------------------------------------------
    for eta in geometric_grid(1e-6, 1e12, 400) {
        let (centered, ratio, d1) = special::bessel_i0_centered_terms(eta);
        emit("bessel_centered_log", eta, centered);
        emit("bessel_ratio", eta, ratio);
        emit("bessel_d1", eta, d1);
        emit(
            "bessel_d2",
            eta,
            special::bessel_i0_centered_second_log_derivative_from_log_abs(eta.ln()),
        );
    }
    // Dense sweep across the ascending/asymptotic seam at 20.
    for eta in linear_grid(15.0, 25.0, 201) {
        let (centered, ratio, d1) = special::bessel_i0_centered_terms(eta);
        emit("bessel_centered_log", eta, centered);
        emit("bessel_ratio", eta, ratio);
        emit("bessel_d1", eta, d1);
        emit(
            "bessel_d2",
            eta,
            special::bessel_i0_centered_second_log_derivative_from_log_abs(eta.ln()),
        );
    }

    // ---- Polygamma family ------------------------------------------------
    for x in geometric_grid(1e-8, 1e10, 400) {
        emit("digamma", x, special::digamma(x));
        emit("trigamma", x, special::trigamma(x));
        emit("tetragamma", x, special::tetragamma(x));
        emit("pentagamma", x, special::pentagamma(x));
    }
    for x in linear_grid(0.5, 30.0, 300) {
        emit("digamma", x, special::digamma(x));
        emit("trigamma", x, special::trigamma(x));
        emit("tetragamma", x, special::tetragamma(x));
        emit("pentagamma", x, special::pentagamma(x));
    }

    // ---- Normal-distribution channels ------------------------------------
    let mut normal_args: Vec<f64> = linear_grid(-40.0, 40.0, 801);
    normal_args.extend(geometric_grid(1e-8, 1e2, 200));
    normal_args.extend(geometric_grid(1e-8, 1e2, 200).into_iter().map(|v| -v));
    for x in normal_args {
        emit("normal_pdf", x, prob::normal_pdf(x));
        emit("normal_cdf", x, prob::normal_cdf(x));
        // The upper tail as its own quantity. `1 - normal_cdf(x)` is exactly
        // zero above x ~ 8.3 and 7% high already at x = 8 (#2562); this channel
        // exists so that cannot silently return.
        emit("normal_sf", x, prob::normal_sf(x));
        emit("normal_logcdf", x, prob::normal_logcdf(x));
        emit("normal_logsf", x, prob::normal_logsf(x));
        let (log_cdf, mills) = prob::signed_probit_logcdf_and_mills_ratio(x);
        emit("probit_logcdf", x, log_cdf);
        emit("probit_mills", x, mills);
        let derivatives = prob::normal_logcdf_derivatives(x);
        for (order, value) in derivatives.iter().enumerate() {
            emit(&format!("logcdf_d{order}"), x, *value);
        }
        if x >= 0.0 {
            emit("erfcx", x, prob::erfcx_nonnegative(x));
        }
    }
    for x in geometric_grid(1e-3, 1e6, 400) {
        emit("erfcx", x, prob::erfcx_nonnegative(x));
    }

    // ---- log1mexp --------------------------------------------------------
    for a in geometric_grid(1e-14, 50.0, 300) {
        emit("log1mexp", a, prob::log1mexp_positive(a));
    }

    // ---- Normal quantile -------------------------------------------------
    for p in geometric_grid(1e-300, 0.5, 400) {
        if let Ok(q) = prob::standard_normal_quantile(p) {
            emit("normal_quantile", p, q);
        }
    }
    for p in linear_grid(0.001, 0.999, 400) {
        if let Ok(q) = prob::standard_normal_quantile(p) {
            emit("normal_quantile", p, q);
        }
    }
    for log_p in linear_grid(-700.0, -1e-6, 400) {
        if let Ok(q) = prob::standard_normal_quantile_from_log_cdf(log_p) {
            emit("normal_quantile_from_log", log_p, q);
        }
    }

    // ---- Gauss-Legendre --------------------------------------------------
    for n in [
        3usize, 4, 5, 7, 8, 12, 15, 16, 20, 24, 31, 32, 40, 48, 63, 64, 80, 96, 100, 127, 128, 160,
        200, 256,
    ] {
        let (nodes, weights) = special::gauss_legendre(n);
        for (i, (node, weight)) in nodes.iter().zip(weights.iter()).enumerate() {
            println!("gl_node\t{n}\t{i}\t{node:e}");
            println!("gl_weight\t{n}\t{i}\t{weight:e}");
        }
    }

    // ---- Binomial coefficient --------------------------------------------
    for n in 0usize..=60 {
        for k in 0..=n {
            println!(
                "binomial\t{n}\t{k}\t{:e}",
                special::binomial_coefficient_f64(n, k)
            );
        }
    }

    // Beta quantiles. The shapes are the ones a beta-regression predictive
    // interval produces from a mean and a variance, plus a direct sweep of
    // small and large shape pairs. The lower tail here reaches quantiles far
    // below `f64::EPSILON`, which is exactly where a solver with an absolute
    // tolerance in `x` stalls rather than degrading (#2528), so the sweep is
    // deliberately weighted toward small `a`.
    for (mu, variance_fraction) in [
        (0.001_f64, 0.3_f64),
        (0.01, 0.2),
        (0.01, 0.5),
        (0.02, 0.3),
        (0.05, 0.5),
        (0.1, 0.5),
        (0.3, 0.5),
        (0.5, 0.5),
        (0.7, 0.5),
        (0.9, 0.2),
    ] {
        let bernoulli_variance = mu * (1.0 - mu);
        let precision = 1.0 / variance_fraction - 1.0;
        let (a, b) = (mu * precision, (1.0 - mu) * precision);
        for p in [0.001_f64, 0.025, 0.1, 0.5, 0.9, 0.975, 0.999] {
            println!(
                "beta_quantile\t{a:e}\t{b:e}\t{p:e}\t{:e}",
                prob::beta_quantile(p, a, b)
            );
        }
        // Keep the moment-matched variance in the record so a reader can see
        // which mean produced which shape pair.
        println!("beta_shape\t{mu:e}\t{bernoulli_variance:e}\t{a:e}\t{b:e}");
    }
    for (a, b) in [
        (0.1_f64, 0.1_f64),
        (0.5, 0.5),
        (0.5, 20.0),
        (20.0, 0.5),
        (1.0, 1.0),
        (2.0, 3.0),
        (2.5, 7.5),
        (100.0, 100.0),
        (1000.0, 5.0),
        (5.0, 1000.0),
    ] {
        for p in [
            1.0e-8_f64, 1.0e-4, 0.001, 0.01, 0.025, 0.1, 0.25, 0.5, 0.75, 0.9, 0.975, 0.999,
        ] {
            println!(
                "beta_quantile\t{a:e}\t{b:e}\t{p:e}\t{:e}",
                prob::beta_quantile(p, a, b)
            );
        }
    }

    // ---- Student-t survival function -------------------------------------
    // Two-argument channel (nu, t). The complement `1 - cdf` saturates SOONER
    // the larger nu is -- already at nu = 500, t = 10, where the true tail is
    // 6.9e-22 -- because the loss is in the subtraction, not in the cdf (#2562).
    for nu in [1.0_f64, 2.5, 5.0, 30.0, 120.0, 500.0, 5000.0, 1.0e4] {
        for t in [
            0.0_f64, 0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 30.0, 40.0, 80.0,
        ] {
            println!(
                "students_t_sf\t{nu:e}\t{t:e}\t{:e}",
                prob::student_t_sf(t, nu)
            );
            println!(
                "students_t_sf\t{nu:e}\t{:e}\t{:e}",
                -t,
                prob::student_t_sf(-t, nu)
            );
        }
    }
}