pub fn normal_cdf(x: f64) -> f64Expand description
Standard normal CDF Phi(x) evaluated via the exact special-function identity
Phi(x) = 0.5 * erfc(-x / sqrt(2)).
This is the exact Gaussian CDF semantics used throughout the codebase. The
numerical erfc implementation may use internal approximations, but the
returned function is the standard normal CDF itself rather than a separate
polynomial surrogate surface.
Examples found in repository?
examples/special_audit_dump.rs (line 79)
33fn main() {
34 // ---- Bessel channels -------------------------------------------------
35 for eta in geometric_grid(1e-6, 1e12, 400) {
36 let (centered, ratio, d1) = special::bessel_i0_centered_terms(eta);
37 emit("bessel_centered_log", eta, centered);
38 emit("bessel_ratio", eta, ratio);
39 emit("bessel_d1", eta, d1);
40 emit(
41 "bessel_d2",
42 eta,
43 special::bessel_i0_centered_second_log_derivative_from_log_abs(eta.ln()),
44 );
45 }
46 // Dense sweep across the ascending/asymptotic seam at 20.
47 for eta in linear_grid(15.0, 25.0, 201) {
48 let (centered, ratio, d1) = special::bessel_i0_centered_terms(eta);
49 emit("bessel_centered_log", eta, centered);
50 emit("bessel_ratio", eta, ratio);
51 emit("bessel_d1", eta, d1);
52 emit(
53 "bessel_d2",
54 eta,
55 special::bessel_i0_centered_second_log_derivative_from_log_abs(eta.ln()),
56 );
57 }
58
59 // ---- Polygamma family ------------------------------------------------
60 for x in geometric_grid(1e-8, 1e10, 400) {
61 emit("digamma", x, special::digamma(x));
62 emit("trigamma", x, special::trigamma(x));
63 emit("tetragamma", x, special::tetragamma(x));
64 emit("pentagamma", x, special::pentagamma(x));
65 }
66 for x in linear_grid(0.5, 30.0, 300) {
67 emit("digamma", x, special::digamma(x));
68 emit("trigamma", x, special::trigamma(x));
69 emit("tetragamma", x, special::tetragamma(x));
70 emit("pentagamma", x, special::pentagamma(x));
71 }
72
73 // ---- Normal-distribution channels ------------------------------------
74 let mut normal_args: Vec<f64> = linear_grid(-40.0, 40.0, 801);
75 normal_args.extend(geometric_grid(1e-8, 1e2, 200));
76 normal_args.extend(geometric_grid(1e-8, 1e2, 200).into_iter().map(|v| -v));
77 for x in normal_args {
78 emit("normal_pdf", x, prob::normal_pdf(x));
79 emit("normal_cdf", x, prob::normal_cdf(x));
80 emit("normal_logcdf", x, prob::normal_logcdf(x));
81 emit("normal_logsf", x, prob::normal_logsf(x));
82 let (log_cdf, mills) = prob::signed_probit_logcdf_and_mills_ratio(x);
83 emit("probit_logcdf", x, log_cdf);
84 emit("probit_mills", x, mills);
85 let derivatives = prob::normal_logcdf_derivatives(x);
86 for (order, value) in derivatives.iter().enumerate() {
87 emit(&format!("logcdf_d{order}"), x, *value);
88 }
89 if x >= 0.0 {
90 emit("erfcx", x, prob::erfcx_nonnegative(x));
91 }
92 }
93 for x in geometric_grid(1e-3, 1e6, 400) {
94 emit("erfcx", x, prob::erfcx_nonnegative(x));
95 }
96
97 // ---- log1mexp --------------------------------------------------------
98 for a in geometric_grid(1e-14, 50.0, 300) {
99 emit("log1mexp", a, prob::log1mexp_positive(a));
100 }
101
102 // ---- Normal quantile -------------------------------------------------
103 for p in geometric_grid(1e-300, 0.5, 400) {
104 if let Ok(q) = prob::standard_normal_quantile(p) {
105 emit("normal_quantile", p, q);
106 }
107 }
108 for p in linear_grid(0.001, 0.999, 400) {
109 if let Ok(q) = prob::standard_normal_quantile(p) {
110 emit("normal_quantile", p, q);
111 }
112 }
113 for log_p in linear_grid(-700.0, -1e-6, 400) {
114 if let Ok(q) = prob::standard_normal_quantile_from_log_cdf(log_p) {
115 emit("normal_quantile_from_log", log_p, q);
116 }
117 }
118
119 // ---- Gauss-Legendre --------------------------------------------------
120 for n in [
121 3usize, 4, 5, 7, 8, 12, 15, 16, 20, 24, 31, 32, 40, 48, 63, 64, 80, 96, 100, 127, 128, 160,
122 200, 256,
123 ] {
124 let (nodes, weights) = special::gauss_legendre(n);
125 for (i, (node, weight)) in nodes.iter().zip(weights.iter()).enumerate() {
126 println!("gl_node\t{n}\t{i}\t{node:e}");
127 println!("gl_weight\t{n}\t{i}\t{weight:e}");
128 }
129 }
130
131 // ---- Binomial coefficient --------------------------------------------
132 for n in 0usize..=60 {
133 for k in 0..=n {
134 println!(
135 "binomial\t{n}\t{k}\t{:e}",
136 special::binomial_coefficient_f64(n, k)
137 );
138 }
139 }
140}