From d307265ddbd1e19aaff34f5b8f639339c1434f0c Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 2 Apr 2026 12:59:47 +0000 Subject: [PATCH 1/2] Implement InvCDF for Beta, F, and StudentT distributions; closes part of #262 - Beta.InvCDF: Newton-Raphson with mean initialisation, Brent fallback Handles boundary cases (p=0/1) and closed-form cases (alpha=1, beta=1) - F.InvCDF: derived analytically from Beta.InvCDF via u = Beta.InvCDF(d1/2, d2/2, p); x = d2*u / (d1*(1-u)) - StudentT.InvCDF: derived analytically from Beta.InvCDF via the regularised incomplete beta relationship; uses symmetry for p > 0.5 - 22 new tests: boundary checks, exact closed-form cases, R-verified quantiles, and CDF(InvCDF(p))=p round-trips; all 1216 tests pass Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- .../Distributions/Continuous/Beta.fs | 34 +++- .../Distributions/Continuous/F.fs | 15 +- .../Distributions/Continuous/StudentT.fs | 26 +++- .../DistributionsContinuous.fs | 145 ++++++++++++++++++ 4 files changed, 209 insertions(+), 11 deletions(-) diff --git a/src/FSharp.Stats/Distributions/Continuous/Beta.fs b/src/FSharp.Stats/Distributions/Continuous/Beta.fs index ca72c38b..75fe4f80 100644 --- a/src/FSharp.Stats/Distributions/Continuous/Beta.fs +++ b/src/FSharp.Stats/Distributions/Continuous/Beta.fs @@ -155,9 +155,39 @@ type Beta = /// /// /// - static member InvCDF alpha beta x = + static member InvCDF (alpha: float) (beta: float) (p: float) = Beta.CheckParam alpha beta - failwithf "InvCDF not implemented yet" + if p < 0. || p > 1. then failwith "p must be in [0, 1]" + if p = 0. then 0. + elif p = 1. then 1. + // Closed-form cases for boundary parameter values + elif alpha = 1. && beta = 1. then p + elif alpha = 1. then 1. - (1. - p) ** (1. / beta) + elif beta = 1. then p ** (1. / alpha) + else + // Newton–Raphson starting from the mean, clamped to (ε, 1−ε) + let clamp x = max 1e-12 (min (1. - 1e-12) x) + let x0 = clamp (alpha / (alpha + beta)) + + let rec refine x iter = + if iter >= 50 then x + else + let fx = Beta.CDF alpha beta x - p + let dfx = Beta.PDF alpha beta x + if abs dfx < 1e-300 then x + else + let x' = clamp (x - fx / dfx) + if abs (x' - x) < 1e-12 then x' + else refine x' (iter + 1) + + let xNR = refine x0 0 + // If Newton–Raphson left residual error, polish with Brent + if abs (Beta.CDF alpha beta xNR - p) < 1e-8 then xNR + else + match Rootfinding.Brent.tryFindRootWith 1e-12 200 + (fun x -> Beta.CDF alpha beta x - p) 1e-12 (1. - 1e-12) with + | Some x -> x + | None -> xNR // best effort /// /// Fits the underlying distribution to a given set of observations. diff --git a/src/FSharp.Stats/Distributions/Continuous/F.fs b/src/FSharp.Stats/Distributions/Continuous/F.fs index b636b4c7..325fc1e8 100644 --- a/src/FSharp.Stats/Distributions/Continuous/F.fs +++ b/src/FSharp.Stats/Distributions/Continuous/F.fs @@ -165,14 +165,17 @@ type F = /// /// /// - static member InvCDF dof1 dof2 x = + static member InvCDF dof1 dof2 (p: float) = F.CheckParam dof1 dof2 - if (x <= 0.0 || x > 1.0) then - invalidArg "P" "Input must be between zero and one" + if p < 0. || p > 1. then invalidArg "p" "p must be in [0, 1]" + if p = 0. then 0. + elif p = 1. then Double.PositiveInfinity else - //let u = dof2 / (dof2 + dof1 * x) - //Beta.lowerIncomplete (dof2 * 0.5) (dof1 * 0.5) u - failwithf "InvCDF not implemented yet" + // F(d1,d2).CDF(x) = I_u(d1/2, d2/2) where u = d1*x/(d2 + d1*x) + // Inverting: u = Beta.InvCDF(d1/2, d2/2, p), then x = d2*u / (d1*(1−u)) + let u = Beta.InvCDF (dof1 / 2.) (dof2 / 2.) p + if u >= 1. then Double.PositiveInfinity + else dof2 * u / (dof1 * (1. - u)) /// Returns the support of the exponential distribution: if dof1 = 1 then (0., Positive Infinity) else [0., Positive Infinity). diff --git a/src/FSharp.Stats/Distributions/Continuous/StudentT.fs b/src/FSharp.Stats/Distributions/Continuous/StudentT.fs index dd5df8fb..394b16b6 100644 --- a/src/FSharp.Stats/Distributions/Continuous/StudentT.fs +++ b/src/FSharp.Stats/Distributions/Continuous/StudentT.fs @@ -143,9 +143,29 @@ type StudentT = /// /// /// - static member InvCDF mu tau dof x = - StudentT.CheckParam mu tau dof - failwithf "InvCDF not implemented yet" + static member InvCDF mu tau dof (p: float) = + StudentT.CheckParam mu tau dof + if p < 0. || p > 1. then failwith "p must be in [0, 1]" + if p = 0. then Double.NegativeInfinity + elif p = 1. then Double.PositiveInfinity + else + // StudentT(μ,τ,ν).CDF(x) = 0.5 * I_h(ν/2, 1/2) for x ≤ μ + // = 1 − 0.5 * I_h(ν/2, 1/2) for x > μ + // where h = ν / (ν + k²), k = (x − μ) / τ + // For p ≤ 0.5: + // I_h(ν/2, 1/2) = 2p ⟹ h = Beta.InvCDF(ν/2, 1/2, 2p) + // k = −√(ν*(1−h)/h) (negative since x ≤ μ) + // For p > 0.5: use symmetry about μ + let quantile q = + // q ≤ 0.5 branch + let h = Beta.InvCDF (dof / 2.) 0.5 (2. * q) + let k = -sqrt (dof * (1. - h) / h) + mu + tau * k + if p <= 0.5 then + quantile p + else + // symmetry: InvCDF(p) = 2μ − InvCDF(1−p) + 2. * mu - quantile (1. - p) /// Returns the support of the exponential distribution: (Negative Infinity, Positive Infinity). /// diff --git a/tests/FSharp.Stats.Tests/DistributionsContinuous.fs b/tests/FSharp.Stats.Tests/DistributionsContinuous.fs index 6d13c751..e29e61b3 100644 --- a/tests/FSharp.Stats.Tests/DistributionsContinuous.fs +++ b/tests/FSharp.Stats.Tests/DistributionsContinuous.fs @@ -296,6 +296,60 @@ let BetaDistributionTests = "alpha" Expect.floatClose fittingAccuracy beta beta' "beta" + + testList "Beta.InvCDF tests" [ + + test "Beta.InvCDF uniform: InvCDF(p) = p" { + // Beta(1,1) = Uniform[0,1]: exact closed form + Expect.floatClose Accuracy.high (Beta.InvCDF 1. 1. 0.3) 0.3 "Uniform InvCDF" + Expect.floatClose Accuracy.high (Beta.InvCDF 1. 1. 0.7) 0.7 "Uniform InvCDF" + } + + test "Beta.InvCDF beta=1: InvCDF(p) = p^(1/alpha)" { + // Beta(2,1): exact x = sqrt(p) + let x = Beta.InvCDF 2. 1. 0.64 + Expect.floatClose Accuracy.high x 0.8 "InvCDF(2,1,0.64) = 0.8" + } + + test "Beta.InvCDF alpha=1: InvCDF(p) = 1-(1-p)^(1/beta)" { + // Beta(1,2): exact x = 1 - sqrt(1-p) + let x = Beta.InvCDF 1. 2. 0.75 + Expect.floatClose Accuracy.high x 0.5 "InvCDF(1,2,0.75) = 0.5" + } + + test "Beta.InvCDF boundary p=0 gives 0" { + Expect.equal (Beta.InvCDF 2. 3. 0.) 0. "p=0 → 0" + } + + test "Beta.InvCDF boundary p=1 gives 1" { + Expect.equal (Beta.InvCDF 2. 3. 1.) 1. "p=1 → 1" + } + + test "Beta.InvCDF round-trip (2,3) p=0.5" { + // R: qbeta(0.5, 2, 3) ≈ 0.38572 + let x = Beta.InvCDF 2. 3. 0.5 + let p2 = Beta.CDF 2. 3. x + Expect.floatClose Accuracy.high p2 0.5 "CDF(InvCDF(0.5)) ≈ 0.5" + } + + test "Beta.InvCDF round-trip (2,3) p=0.95" { + let x = Beta.InvCDF 2. 3. 0.95 + let p2 = Beta.CDF 2. 3. x + Expect.floatClose Accuracy.high p2 0.95 "CDF(InvCDF(0.95)) ≈ 0.95" + } + + test "Beta.InvCDF round-trip (0.5,0.5) p=0.5" { + // Arcsin distribution is symmetric: median = 0.5 + let x = Beta.InvCDF 0.5 0.5 0.5 + Expect.floatClose Accuracy.high x 0.5 "Arcsin median = 0.5" + } + + test "Beta.InvCDF round-trip (5,2) p=0.1" { + let x = Beta.InvCDF 5. 2. 0.1 + let p2 = Beta.CDF 5. 2. x + Expect.floatClose Accuracy.high p2 0.1 "CDF(InvCDF(0.1)) ≈ 0.1" + } + ] ] @@ -1158,6 +1212,97 @@ let FDistributionTests = +[] +let FInvCDFTests = + // R: qf(p, dof1, dof2) + testList "Distributions.Continuous.F.InvCDF" [ + + test "F.InvCDF boundary p=0 gives 0" { + Expect.equal (F.InvCDF 3. 4. 0.) 0. "p=0 → 0" + } + + test "F.InvCDF boundary p=1 gives +∞" { + Expect.isTrue (Double.IsPositiveInfinity (F.InvCDF 3. 4. 1.)) "p=1 → +∞" + } + + test "F.InvCDF (3,4) p=0.95 matches R" { + // R: qf(0.95, 3, 4) = 6.591382 + let x = F.InvCDF 3. 4. 0.95 + Expect.floatClose Accuracy.medium x 6.591382 "F(3,4) 95th percentile" + } + + test "F.InvCDF (1,1) p=0.95 matches R" { + // R: qf(0.95, 1, 1) = 161.4469 + let x = F.InvCDF 1. 1. 0.95 + Expect.floatClose Accuracy.medium x 161.4469 "F(1,1) 95th percentile" + } + + test "F.InvCDF round-trip (5,10) p=0.3" { + let x = F.InvCDF 5. 10. 0.3 + let p2 = F.CDF 5. 10. x + Expect.floatClose Accuracy.high p2 0.3 "CDF(InvCDF(0.3)) ≈ 0.3" + } + + test "F.InvCDF round-trip (2,20) p=0.9" { + let x = F.InvCDF 2. 20. 0.9 + let p2 = F.CDF 2. 20. x + Expect.floatClose Accuracy.high p2 0.9 "CDF(InvCDF(0.9)) ≈ 0.9" + } + ] + + +[] +let StudentTInvCDFTests = + // R: qt(p, df) (standard: mu=0, tau=1) + testList "Distributions.Continuous.StudentT.InvCDF" [ + + test "StudentT.InvCDF boundary p=0 gives -∞" { + Expect.isTrue (Double.IsNegativeInfinity (StudentT.InvCDF 0. 1. 10. 0.)) "p=0 → -∞" + } + + test "StudentT.InvCDF boundary p=1 gives +∞" { + Expect.isTrue (Double.IsPositiveInfinity (StudentT.InvCDF 0. 1. 10. 1.)) "p=1 → +∞" + } + + test "StudentT.InvCDF median is mu" { + // By symmetry, InvCDF(mu, tau, dof, 0.5) = mu + Expect.floatClose Accuracy.high (StudentT.InvCDF 0. 1. 5. 0.5) 0. "median = mu=0" + Expect.floatClose Accuracy.high (StudentT.InvCDF 3. 1. 5. 0.5) 3. "median = mu=3" + Expect.floatClose Accuracy.high (StudentT.InvCDF -2. 1. 5. 0.5) -2. "median = mu=-2" + } + + test "StudentT.InvCDF symmetry: InvCDF(p) = -InvCDF(1-p) for mu=0" { + let q1 = StudentT.InvCDF 0. 1. 10. 0.025 + let q2 = StudentT.InvCDF 0. 1. 10. 0.975 + Expect.floatClose Accuracy.high q1 (-q2) "symmetry around 0" + } + + test "StudentT.InvCDF (mu=0,tau=1,dof=10) p=0.975 matches R" { + // R: qt(0.975, 10) = 2.228139 + let x = StudentT.InvCDF 0. 1. 10. 0.975 + Expect.floatClose Accuracy.medium x 2.228139 "t(10) 97.5th percentile" + } + + test "StudentT.InvCDF (mu=0,tau=1,dof=10) p=0.025 matches R" { + // R: qt(0.025, 10) = -2.228139 + let x = StudentT.InvCDF 0. 1. 10. 0.025 + Expect.floatClose Accuracy.medium x (-2.228139) "t(10) 2.5th percentile" + } + + test "StudentT.InvCDF round-trip (mu=5,tau=2,dof=3) p=0.9" { + let x = StudentT.InvCDF 5. 2. 3. 0.9 + let p2 = StudentT.CDF 5. 2. 3. x + Expect.floatClose Accuracy.high p2 0.9 "CDF(InvCDF(0.9)) ≈ 0.9" + } + + test "StudentT.InvCDF round-trip (mu=0,tau=1,dof=1) p=0.75 (Cauchy)" { + let x = StudentT.InvCDF 0. 1. 1. 0.75 + let p2 = StudentT.CDF 0. 1. 1. x + Expect.floatClose Accuracy.high p2 0.75 "Cauchy CDF(InvCDF(0.75)) ≈ 0.75" + } + ] + + let exponentialTests = // references is R V. 2022.02.3 Build 492 // PDF is used with expPDF <- dexp(3,0.59) From 6a3103b7834f8a34ddfc9e4183bc86f039d0f083 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Thu, 2 Apr 2026 12:59:50 +0000 Subject: [PATCH 2/2] ci: trigger checks