spb/metrika Public

1//2//  Distributions.swift3//  Metrika4//5//  Author:  Simon-Pierre Boucher6//  Contact: contact@spboucher.ai7//  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.8//910import Foundation1112/// Special functions and distribution CDFs used for inference output.13/// Implementations follow the classic Lanczos / continued-fraction14/// formulations and are accurate to well beyond the 1e-10 project15/// tolerance for the parameter ranges used in regression output.16public enum ZQDistributions {1718    /// Natural log of the gamma function (Lanczos approximation, g=7).19    public static func logGamma(_ x: Double) -> Double {20        precondition(x > 0, "logGamma requires x > 0")21        let coefficients: [Double] = [22            676.5203681218851, -1259.1392167224028, 771.32342877765313,23            -176.61502916214059, 12.507343278686905, -0.13857109526572012,24            9.9843695780195716e-6, 1.5056327351493116e-7,25        ]26        if x < 0.5 {27            // Reflection formula.28            return log(.pi / sin(.pi * x)) - logGamma(1 - x)29        }30        let z = x - 131        var sum = 0.9999999999998099332        for (i, c) in coefficients.enumerated() {33            sum += c / (z + Double(i) + 1)34        }35        let t = z + 7.536        return 0.5 * log(2 * .pi) + (z + 0.5) * log(t) - t + log(sum)37    }3839    /// Regularized incomplete beta function I_x(a, b) via the Lentz40    /// continued-fraction algorithm.41    public static func incompleteBeta(a: Double, b: Double, x: Double) -> Double {42        precondition(a > 0 && b > 0, "incompleteBeta requires a, b > 0")43        if x <= 0 { return 0 }44        if x >= 1 { return 1 }4546        let logBeta = logGamma(a + b) - logGamma(a) - logGamma(b)47        let front = exp(logBeta + a * log(x) + b * log(1 - x))4849        // Use the symmetry relation to keep the continued fraction in its50        // rapidly converging region.51        if x < (a + 1) / (a + b + 2) {52            return front * betaContinuedFraction(a: a, b: b, x: x) / a53        } else {54            return 1 - front * betaContinuedFraction(a: b, b: a, x: 1 - x) / b55        }56    }5758    private static func betaContinuedFraction(a: Double, b: Double, x: Double) -> Double {59        let tiny = 1e-30060        let epsilon = 1e-1661        let qab = a + b62        let qap = a + 163        let qam = a - 16465        var c = 1.066        var d = 1 - qab * x / qap67        if abs(d) < tiny { d = tiny }68        d = 1 / d69        var h = d7071        for m in 1...300 {72            let dm = Double(m)73            // Even step.74            var numerator = dm * (b - dm) * x / ((qam + 2 * dm) * (a + 2 * dm))75            d = 1 + numerator * d76            if abs(d) < tiny { d = tiny }77            c = 1 + numerator / c78            if abs(c) < tiny { c = tiny }79            d = 1 / d80            h *= d * c81            // Odd step.82            numerator = -(a + dm) * (qab + dm) * x / ((a + 2 * dm) * (qap + 2 * dm))83            d = 1 + numerator * d84            if abs(d) < tiny { d = tiny }85            c = 1 + numerator / c86            if abs(c) < tiny { c = tiny }87            d = 1 / d88            let delta = d * c89            h *= delta90            if abs(delta - 1) < epsilon { break }91        }92        return h93    }9495    /// Standard normal CDF.96    public static func normalCDF(_ z: Double) -> Double {97        0.5 * erfc(-z / 2.0.squareRoot())98    }99100    /// Student t CDF with `df` degrees of freedom.101    public static func studentTCDF(_ t: Double, df: Double) -> Double {102        precondition(df > 0, "studentTCDF requires df > 0")103        let x = df / (df + t * t)104        let tail = 0.5 * incompleteBeta(a: df / 2, b: 0.5, x: x)105        return t > 0 ? 1 - tail : tail106    }107108    /// Two-sided p-value for a t statistic, computed directly from the109    /// tail integral: 2·P(T > |t|) = I_{df/(df+t²)}(df/2, 1/2). Going110    /// through `2·(1 − CDF)` would lose all precision below ~1e-16.111    public static func tTestPValue(_ t: Double, df: Double) -> Double {112        incompleteBeta(a: df / 2, b: 0.5, x: df / (df + t * t))113    }114115    /// F distribution CDF with (df1, df2) degrees of freedom.116    public static func fCDF(_ f: Double, df1: Double, df2: Double) -> Double {117        precondition(df1 > 0 && df2 > 0, "fCDF requires positive df")118        if f <= 0 { return 0 }119        let x = df1 * f / (df1 * f + df2)120        return incompleteBeta(a: df1 / 2, b: df2 / 2, x: x)121    }122123    /// Upper-tail p-value for an F statistic, computed directly as124    /// I_{df2/(df2+df1·f)}(df2/2, df1/2) to preserve precision for large F.125    public static func fTestPValue(_ f: Double, df1: Double, df2: Double) -> Double {126        if f <= 0 { return 1 }127        return incompleteBeta(a: df2 / 2, b: df1 / 2, x: df2 / (df2 + df1 * f))128    }129130    /// Regularized lower incomplete gamma P(a, x), via the series for131    /// x < a+1 and the Lentz continued fraction otherwise.132    public static func incompleteGamma(a: Double, x: Double) -> Double {133        precondition(a > 0, "incompleteGamma requires a > 0")134        if x <= 0 { return 0 }135136        if x < a + 1 {137            // Series representation.138            var term = 1.0 / a139            var sum = term140            var ap = a141            for _ in 0..<500 {142                ap += 1143                term *= x / ap144                sum += term145                if abs(term) < abs(sum) * 1e-16 { break }146            }147            return sum * exp(-x + a * log(x) - logGamma(a))148        } else {149            // Continued fraction for Q(a, x); P = 1 − Q.150            let tiny = 1e-300151            var b = x + 1 - a152            var c = 1 / tiny153            var d = 1 / b154            var h = d155            for i in 1...500 {156                let an = -Double(i) * (Double(i) - a)157                b += 2158                d = an * d + b159                if abs(d) < tiny { d = tiny }160                c = b + an / c161                if abs(c) < tiny { c = tiny }162                d = 1 / d163                let delta = d * c164                h *= delta165                if abs(delta - 1) < 1e-16 { break }166            }167            let q = exp(-x + a * log(x) - logGamma(a)) * h168            return 1 - q169        }170    }171172    /// Chi-square CDF with `df` degrees of freedom.173    public static func chiSquareCDF(_ x: Double, df: Double) -> Double {174        x <= 0 ? 0 : incompleteGamma(a: df / 2, x: x / 2)175    }176177    /// Upper-tail p-value for a chi-square statistic.178    public static func chiSquarePValue(_ x: Double, df: Double) -> Double {179        1 - chiSquareCDF(x, df: df)180    }181182    /// Standard normal quantile via bisection on the CDF (erfc-based, so183    /// accurate to ~1e-15) — used for z confidence intervals.184    public static func normalQuantile(_ p: Double) -> Double {185        precondition(p > 0 && p < 1, "quantile requires 0 < p < 1")186        if abs(p - 0.5) < 1e-15 { return 0 }187        var low = -40.0, high = 40.0188        for _ in 0..<200 {189            let mid = 0.5 * (low + high)190            if normalCDF(mid) < p { low = mid } else { high = mid }191            if high - low < 1e-15 * max(1, abs(mid)) { break }192        }193        return 0.5 * (low + high)194    }195196    /// Student t quantile (inverse CDF) via bisection refined with Newton197    /// steps — used for confidence intervals.198    public static func studentTQuantile(_ p: Double, df: Double) -> Double {199        precondition(p > 0 && p < 1, "quantile requires 0 < p < 1")200        if abs(p - 0.5) < 1e-15 { return 0 }201202        var low = -1000.0, high = 1000.0203        var mid = 0.0204        for _ in 0..<200 {205            mid = 0.5 * (low + high)206            if studentTCDF(mid, df: df) < p {207                low = mid208            } else {209                high = mid210            }211            if high - low < 1e-13 * max(1, abs(mid)) { break }212        }213        return 0.5 * (low + high)214    }215}216