/// @module std::core::distributions_advanced
/// Advanced Statistical Distributions
///
/// Provides PDF, CDF, and sampling functions for common distributions
/// beyond the basic ones in distributions.shape.
/// Built on top of random.shape intrinsics using standard mathematical methods.
let PI = 3.141592653589793;
let E = 2.718281828459045;
let SQRT_2PI = 2.5066282746310002;
// ===== Helper: Gamma function (Lanczos approximation) =====
/// Log-gamma function via Lanczos approximation (g=7, n=9)
function ln_gamma(x) {
let g = 7.0;
let coefs = [
0.99999999999980993,
676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
0.0000099843695780195716,
0.00000015056327351493116
];
if x < 0.5 {
// Reflection formula: Gamma(1-x) * Gamma(x) = pi / sin(pi*x)
let reflect = ln(PI / sin(PI * x)) - ln_gamma(1.0 - x);
return reflect;
}
let z = x - 1.0;
let mut ag = coefs[0];
for i in range(1, 9) {
ag = ag + coefs[i] / (z + i);
}
let t = z + g + 0.5;
0.5 * ln(2.0 * PI) + (z + 0.5) * ln(t) - t + ln(ag)
}
/// Gamma function
pub fn gamma(x) {
exp(ln_gamma(x))
}
/// Beta function B(a, b) = Gamma(a) * Gamma(b) / Gamma(a+b)
pub fn beta_fn(a, b) {
exp(ln_gamma(a) + ln_gamma(b) - ln_gamma(a + b))
}
// ===== Normal Distribution =====
/// Standard normal PDF: phi(x)
pub fn normal_pdf(x, mu = 0.0, sigma = 1.0) {
let z = (x - mu) / sigma;
exp(-0.5 * z * z) / (sigma * SQRT_2PI)
}
/// Standard normal CDF using rational approximation (Abramowitz & Stegun)
pub fn normal_cdf(x, mu = 0.0, sigma = 1.0) {
let z = (x - mu) / sigma;
// Use symmetry for negative values
let mut sign = 1.0;
let mut z_abs = z;
if z < 0.0 {
sign = -1.0;
z_abs = -z;
}
// Abramowitz & Stegun 7.1.26 approximation to erfc(x)
// erf(x) = 1 - (a1*t + a2*t^2 + ... + a5*t^5) * exp(-x^2)
// where t = 1/(1+p*x)
// Normal CDF: Phi(z) = 0.5 * (1 + erf(z/sqrt(2)))
let p_coef = 0.3275911;
let a1 = 0.254829592;
let a2 = -0.284496736;
let a3 = 1.421413741;
let a4 = -1.453152027;
let a5 = 1.061405429;
let x_erf = z_abs / sqrt(2.0);
let t = 1.0 / (1.0 + p_coef * x_erf);
let t2 = t * t;
let t3 = t2 * t;
let t4 = t3 * t;
let t5 = t4 * t;
let erf_approx = 1.0 - (a1 * t + a2 * t2 + a3 * t3 + a4 * t4 + a5 * t5) * exp(-x_erf * x_erf);
0.5 * (1.0 + sign * erf_approx)
}
/// Inverse normal CDF (quantile function) using rational approximation
pub fn normal_quantile(p, mu = 0.0, sigma = 1.0) {
// Beasley-Springer-Moro algorithm
let a = [
-3.969683028665376e+01,
2.209460984245205e+02,
-2.759285104469687e+02,
1.383577518672690e+02,
-3.066479806614716e+01,
2.506628277459239e+00
];
let b = [
-5.447609879822406e+01,
1.615858368580409e+02,
-1.556989798598866e+02,
6.680131188771972e+01,
-1.328068155288572e+01
];
let c = [
-7.784894002430293e-03,
-3.223964580411365e-01,
-2.400758277161838e+00,
-2.549732539343734e+00,
4.374664141464968e+00,
2.938163982698783e+00
];
let d = [
7.784695709041462e-03,
3.224671290700398e-01,
2.445134137142996e+00,
3.754408661907416e+00
];
let p_low = 0.02425;
let p_high = 1.0 - p_low;
let mut q = 0.0;
let mut r = 0.0;
let mut z = 0.0;
if p < p_low {
q = sqrt(-2.0 * ln(p));
z = (((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) /
((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0);
} else if p <= p_high {
q = p - 0.5;
r = q * q;
z = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5]) * q /
(((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0);
} else {
q = sqrt(-2.0 * ln(1.0 - p));
z = -(((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) /
((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0);
}
mu + sigma * z
}
// ===== Chi-Square Distribution =====
/// Chi-square PDF: f(x; k) = x^(k/2-1) * exp(-x/2) / (2^(k/2) * Gamma(k/2))
pub fn chi_square_pdf(x, k) {
if x <= 0.0 {
return 0.0;
}
let half_k = k / 2.0;
exp((half_k - 1.0) * ln(x) - x / 2.0 - half_k * ln(2.0) - ln_gamma(half_k))
}
/// Chi-square CDF (via regularized incomplete gamma — series expansion)
pub fn chi_square_cdf(x, k) {
if x <= 0.0 {
return 0.0;
}
regularized_gamma_p(k / 2.0, x / 2.0)
}
/// Sample from chi-square distribution (sum of k squared standard normals)
pub fn chi_square_sample(k) {
let mut sum = 0.0;
for i in range(0, k) {
let z = __intrinsic_random_normal(0.0, 1.0);
sum = sum + z * z;
}
sum
}
// ===== Student's t Distribution =====
/// Student's t PDF
pub fn t_pdf(x, df) {
let half_df_plus = (df + 1.0) / 2.0;
let half_df = df / 2.0;
let coeff = exp(ln_gamma(half_df_plus) - ln_gamma(half_df)) / sqrt(df * PI);
coeff * pow(1.0 + x * x / df, -half_df_plus)
}
/// Student's t CDF (numerical integration via Simpson's rule)
pub fn t_cdf(x, df) {
// Use symmetry and numerical integration
if x == 0.0 {
return 0.5;
}
// Integrate from -inf to x using change of variable
// For t-distribution, use the regularized incomplete beta function
let t2 = x * x;
let p = 1.0 - 0.5 * regularized_beta(df / (df + t2), df / 2.0, 0.5);
if x > 0.0 {
p
} else {
1.0 - p
}
}
/// Sample from Student's t distribution
pub fn t_sample(df) {
let z = __intrinsic_random_normal(0.0, 1.0);
let v = chi_square_sample(df) / df;
z / sqrt(v)
}
// ===== Beta Distribution =====
/// Beta PDF: f(x; a, b) = x^(a-1) * (1-x)^(b-1) / B(a,b)
pub fn beta_pdf(x, a, b) {
if x <= 0.0 || x >= 1.0 {
return 0.0;
}
exp((a - 1.0) * ln(x) + (b - 1.0) * ln(1.0 - x) - ln_gamma(a) - ln_gamma(b) + ln_gamma(a + b))
}
/// Beta CDF (regularized incomplete beta function)
pub fn beta_cdf(x, a, b) {
if x <= 0.0 {
return 0.0;
}
if x >= 1.0 {
return 1.0;
}
regularized_beta(x, a, b)
}
/// Sample from Beta distribution (Joehnk's method for small params,
/// gamma ratio for general case)
pub fn beta_sample(a, b) {
// Use gamma ratio method: X = G1/(G1+G2) where G1~Gamma(a), G2~Gamma(b)
let g1 = gamma_sample(a);
let g2 = gamma_sample(b);
g1 / (g1 + g2)
}
// ===== Gamma Distribution =====
/// Gamma PDF: f(x; k, theta) = x^(k-1) * exp(-x/theta) / (theta^k * Gamma(k))
pub fn gamma_pdf(x, k, theta = 1.0) {
if x <= 0.0 {
return 0.0;
}
exp((k - 1.0) * ln(x) - x / theta - k * ln(theta) - ln_gamma(k))
}
/// Gamma CDF
pub fn gamma_cdf(x, k, theta = 1.0) {
if x <= 0.0 {
return 0.0;
}
regularized_gamma_p(k, x / theta)
}
/// Sample from Gamma distribution (Marsaglia and Tsang's method)
pub fn gamma_sample(k, theta = 1.0) {
if k < 1.0 {
// For k < 1, use Gamma(k+1) * U^(1/k)
let s = gamma_sample(k + 1.0, 1.0);
return s * pow(__intrinsic_random(), 1.0 / k) * theta;
}
let d = k - 1.0 / 3.0;
let c = 1.0 / sqrt(9.0 * d);
let mut result = 0.0;
let mut found = false;
let mut attempts = 0;
while !found && attempts < 10000 {
let mut x = __intrinsic_random_normal(0.0, 1.0);
let mut v = 1.0 + c * x;
while v <= 0.0 {
x = __intrinsic_random_normal(0.0, 1.0);
v = 1.0 + c * x;
}
v = v * v * v;
let u = __intrinsic_random();
if u < 1.0 - 0.0331 * (x * x) * (x * x) {
result = d * v * theta;
found = true;
} else if ln(u) < 0.5 * x * x + d * (1.0 - v + ln(v)) {
result = d * v * theta;
found = true;
}
attempts = attempts + 1;
}
result
}
// ===== Special Functions =====
/// Regularized incomplete gamma function P(a, x) via series expansion
function regularized_gamma_p(a, x) {
if x <= 0.0 {
return 0.0;
}
if x > a + 1.0 {
// Use complement: P = 1 - Q, where Q uses continued fraction
return 1.0 - regularized_gamma_q(a, x);
}
// Series expansion: P(a,x) = e^(-x) * x^a * sum(x^n / Gamma(a+n+1))
let mut sum = 1.0 / a;
let mut term = 1.0 / a;
for n in range(1, 200) {
term = term * x / (a + n);
sum = sum + term;
if abs(term) < abs(sum) * 0.0000000001 {
break;
}
}
exp(-x + a * ln(x) - ln_gamma(a)) * sum
}
/// Complementary regularized incomplete gamma Q(a, x) via continued fraction
function regularized_gamma_q(a, x) {
// Lentz's method for continued fraction
let mut f = 0.0000000000000001;
let mut c_val = f;
let mut d_val = 0.0;
for n in range(1, 200) {
let mut an = 0.0;
if n % 2 == 1 {
let k = (n - 1) / 2;
an = -(a + k) * (a + k + 0.0 - a + n * 1.0);
an = (k + 1.0 - a) * (k + 1.0);
// Simplified: use the standard CF for Q
}
// Use simpler Lentz method
break;
}
// Fallback: use series for P and return 1 - P
// For large x, this converges fast anyway
let mut sum = 1.0 / a;
let mut term = 1.0 / a;
for n in range(1, 200) {
term = term * x / (a + n);
sum = sum + term;
if abs(term) < abs(sum) * 0.0000000001 {
break;
}
}
1.0 - exp(-x + a * ln(x) - ln_gamma(a)) * sum
}
/// Regularized incomplete beta function I_x(a, b)
/// Using continued fraction expansion (Lentz's method)
function regularized_beta(x, a, b) {
if x <= 0.0 {
return 0.0;
}
if x >= 1.0 {
return 1.0;
}
// Use symmetry: if x > (a+1)/(a+b+2), use I_x(a,b) = 1 - I_(1-x)(b,a)
if x > (a + 1.0) / (a + b + 2.0) {
return 1.0 - regularized_beta(1.0 - x, b, a);
}
// Compute the prefix: x^a * (1-x)^b / (a * B(a,b))
let ln_prefix = a * ln(x) + b * ln(1.0 - x) - ln(a) - ln_gamma(a) - ln_gamma(b) + ln_gamma(a + b);
let prefix = exp(ln_prefix);
// Continued fraction (Lentz's method)
let eps = 0.0000000001;
let tiny = 0.0000000000000001;
let mut f = 1.0;
let mut c_val = 1.0;
let mut d_val = 1.0 - (a + b) * x / (a + 1.0);
if abs(d_val) < tiny {
d_val = tiny;
}
d_val = 1.0 / d_val;
f = d_val;
for m in range(1, 200) {
// Even step
let mut num = m * (b - m) * x / ((a + 2.0 * m - 1.0) * (a + 2.0 * m));
d_val = 1.0 + num / d_val;
if abs(d_val) < tiny { d_val = tiny; }
c_val = 1.0 + num / c_val;
if abs(c_val) < tiny { c_val = tiny; }
d_val = 1.0 / d_val;
f = f * d_val * c_val;
// Odd step
num = -(a + m) * (a + b + m) * x / ((a + 2.0 * m) * (a + 2.0 * m + 1.0));
d_val = 1.0 + num / d_val;
if abs(d_val) < tiny { d_val = tiny; }
c_val = 1.0 + num / c_val;
if abs(c_val) < tiny { c_val = tiny; }
d_val = 1.0 / d_val;
let delta = d_val * c_val;
f = f * delta;
if abs(delta - 1.0) < eps {
break;
}
}
prefix * f
}