From 55f8f1e44098dbe65d496208a05ef6df5945ed85 Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Wed, 28 May 2025 15:08:53 +0200 Subject: [PATCH 01/12] Update SpecialFunctions/Gamma.fs to support generic type 'T as implemented in FsMath. --- src/FSharp.Stats/SpecialFunctions/Gamma.fs | 505 +++++++++++++-------- 1 file changed, 326 insertions(+), 179 deletions(-) diff --git a/src/FSharp.Stats/SpecialFunctions/Gamma.fs b/src/FSharp.Stats/SpecialFunctions/Gamma.fs index 189af2aa..0e662773 100644 --- a/src/FSharp.Stats/SpecialFunctions/Gamma.fs +++ b/src/FSharp.Stats/SpecialFunctions/Gamma.fs @@ -1,6 +1,7 @@ namespace FSharp.Stats.SpecialFunctions open System +open FSharp.Stats.GenericMath /// Approximations for the gamma function and related functions. /// @@ -23,15 +24,35 @@ module Gamma = /// The caller is responsible to handle edge cases such as nan, infinity, and -infinity in the input /// /// The function input for approximating Γ(z) - let _gamma z = - let lanczosCoefficients = [76.18009172947146;-86.50532032941677;24.01409824083091;-1.231739572450155;0.1208650973866179e-2;-0.5395239384953e-5] - let rec sumCoefficients acc i coefficients = - match coefficients with - | [] -> acc - | h::t -> sumCoefficients (acc + (h/i)) (i+1.0) t - let gamma = 5.0 - let x = z - 1.0 - Math.Pow(x + gamma + 0.5, x + 0.5) * Math.Exp( -(x + gamma + 0.5) ) * Math.Sqrt( 2.0 * Math.PI ) * sumCoefficients 1.000000000190015 (x + 1.0) lanczosCoefficients + let inline _gamma<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> > (z: 'T) : 'T = + let coeffs : Vector<'T> = + [| + T 76.18009172947146 + T -86.50532032941677 + T 24.01409824083091 + T -1.231739572450155 + T 0.001208650973866179 + T -0.000005395239384953 + |] + + let x = z - one + let g = T 5.0 + let half = T 0.5 + let xg = x + g + half + + let sum = + coeffs + |> Vector.foldi (fun i acc c -> + acc + c / (x + T (float (i + 1)))) (T 1.000000000190015) + + pow xg (x + half) + * exp (-xg) + * sqrt ('T.Pi * T 2.0) + * sum /// /// Computes an approximation of the real value of the log gamma function using the Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -40,16 +61,42 @@ module Gamma = /// The caller is responsible to handle edge cases such as nan, infinity, and -infinity in the input /// /// The function input for approximating ln(Γ(z)) - let _gammaLn z = - let lanczosCoefficients = [76.18009172947146;-86.50532032941677;24.01409824083091;-1.231739572450155;0.1208650973866179e-2;-0.5395239384953e-5] - let rec sumCoefficients acc i coefficients = - match coefficients with - | [] -> acc - | h::t -> sumCoefficients (acc + (h/i)) (i+1.0) t - let gamma = 5.0 - let x = z - 1.0 - let tmp = x + gamma + 0.5 - -(tmp - ((x + 0.5) * log(tmp))) + log(Math.Sqrt( 2.0 * Math.PI ) * sumCoefficients 1.000000000190015 (x + 1.0) lanczosCoefficients) + let inline _gammaLn<'T when 'T :> Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> > (z: 'T) : 'T = + let coeffs = + [| + T 76.18009172947146 + T -86.50532032941677 + T 24.01409824083091 + T -1.231739572450155 + T 0.001208650973866179 + T -0.000005395239384953 + |] + + let pi = 'T.Pi + let x = z - one + let g = T 5.0 + let half = T 0.5 + let xg = x + g + half + + // let sum = + // coeffs + // |> Vector.foldi (fun i acc c -> + // acc + c / (x + T (float (i + 1)))) (T 1.000000000190015) + + // //(x + half) * log xg - xg + log (sqrt ('T.Pi * T<'T> 2.0) * sum) + // (x + half) * log xg - xg + T 0.5 * log (T (2.0 * System.Math.PI)) + log sum + let sum = + coeffs + |> Vector.foldi (fun i acc c -> acc + c / (z + T (float i))) (T 1.000000000190015) + + (x + half) * log xg + - xg + + half * log (T 2.0 * pi) + + log sum /// /// Computes an approximation of the real value of the gamma function using the Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -59,12 +106,18 @@ module Gamma = /// This might be slower than the unchecked version `_gamma` but does not require input sanitation to get expected results for these cases. /// /// The function input for approximating Γ(z) - let gamma z = + let inline gamma<'T when 'T :> Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : equality > (z: 'T) : 'T = //TODO: maybe rename in "complete" for consisteny match z with - | z when (infinity.Equals(z)) -> infinity - | z when ((-infinity).Equals(z)) -> nan - | _ -> _gamma z + | z when 'T.IsPositiveInfinity z -> z + | z_ when 'T.IsNegativeInfinity z -> T nan + | z when 'T.IsNaN z -> z + | z -> _gamma z + /// /// Computes an approximation of the real value of the log gamma function using the Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -74,146 +127,193 @@ module Gamma = /// This might be slower than the unchecked version `_gamma` but does not require input sanitation to get expected results for these cases. /// /// The function input for approximating ln(Γ(z)) - let gammaLn z = + let inline gammaLn<'T when 'T :> Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : equality > (z: 'T) : 'T = + match z with - | z when (infinity.Equals(z)) -> infinity - | z when ((-infinity).Equals(z)) -> nan - | _ -> _gammaLn z + | z when 'T.IsPositiveInfinity z -> z + | z when 'T.IsNegativeInfinity z -> T<'T> System.Double.NaN + | z when 'T.IsNaN z -> z + | z -> _gammaLn z let private FPMIN = 1.0e-30 //(System.Double.MinValue + System.Double.Epsilon) let private EPS = 3.0e-8 //System.Double.Epsilon - let private gser a x = - let gln = gammaLn a - let rec loop sum del ap = - let ap' = ap + 1. + let inline gser<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : comparison > (a: 'T) (x: 'T) : 'T = + + // Relative precision threshold used in convergence tests + let EPS = T<'T> 3.0e-8 + + let gln = _gammaLn a + let rec loop sum (del: 'T) (ap: 'T) = + let ap':'T = ap + 'T.One let del' = del * x / ap' - if (abs del < (abs sum) * EPS) then - sum * exp(-x + a * log(x) - gln) - else - loop (sum+del') (del') (ap') - - loop (1. / a) (1. / a) a + if abs del < abs sum * EPS then + sum * exp (-x + a * log x - gln) + else + loop (sum + del') del' ap' + + loop (T<'T> 1.0 / a) (T<'T> 1.0 / a) a // Page 286 // Returns the incomplete gamma function Q(a,x) evaluated by its continued fraction representation. Also sets ln // .a/ as gln. User should not call directly. - let private gcf a x = - - // let eps = 3.0e-7//System.Double.Epsilon - let gln = gammaLn a - let b = x + 1.0 - a - let c = 1. / FPMIN - let d = 1. / b - let h = d - - let rec loop i b c d h = + let inline gcf<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : comparison > (a: 'T) (x: 'T) : 'T = + + // Smallest positive number to avoid division underflow in continued fractions + let FPMIN = T<'T> 1.0e-30 + // Relative precision threshold used in convergence tests + let EPS = T<'T> 3.0e-8 + + let gln = _gammaLn a + let one = 'T.One + let b0 = x + one - a + let c0 = one / FPMIN + let d0 = one / b0 + let h0 = d0 + + let rec loop (i: 'T) (b: 'T) (c: 'T) (d: 'T) (h: 'T) = let an = -i * (i - a) - let b = b + 2. - let d = + let b = b + T<'T> 2.0 + let d = let tmp = an * d + b if abs tmp < FPMIN then FPMIN else tmp - let c = + let c = let tmp = b + an / c if abs tmp < FPMIN then FPMIN else tmp - let d = 1. / d + let d = one / d let del = d * c let h = h * del - if (abs (del - 1.) <= EPS) then - exp(-x + a * log(x) - gln) * h + if abs (del - one) <= EPS then + exp (-x + a * log x - gln) * h else - loop (i+1.) b c d h - - loop (1.) b c d h + loop (i + one) b c d h + + loop one b0 c0 d0 h0 // Incomplete gamma by quadrature. Returns P(a,x) or Q(a,x), when psig is 1 or 0, respectively. // User should not call directly. - let private gammpapprox a x (psig:bool) = + let inline gammpapprox<'T when 'T :> Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : comparison > (a: 'T) (x: 'T) (psig:bool) : 'T = // Abscissas for Gauss-Legendre quadrature - let gaussLegQuadY = - [| - 0.0021695375159141994; 0.011413521097787704; 0.027972308950302116; 0.051727015600492421; 0.082502225484340941; 0.12007019910960293; - 0.16415283300752470; 0.21442376986779355; 0.27051082840644336; 0.33199876341447887; 0.39843234186401943; 0.46931971407375483; - 0.54413605556657973; 0.62232745288031077; 0.70331500465597174; 0.78649910768313447; 0.87126389619061517; 0.95698180152629142; - |] - // Weights for Gauss-Legendre quadrature - let gaussLegQuadWeights = - [| - 0.0055657196642445571; 0.012915947284065419; 0.020181515297735382; 0.027298621498568734; 0.034213810770299537; 0.040875750923643261; - 0.047235083490265582; 0.053244713977759692; 0.058860144245324798; 0.064039797355015485; 0.068745323835736408; 0.072941885005653087; - 0.076598410645870640; 0.079687828912071670; 0.082187266704339706; 0.084078218979661945; 0.085346685739338721; 0.085983275670394821; - |] - - let ngau = gaussLegQuadWeights.Length //18 + let gaussLegQuadY = [| + 0.0021695375159141994; 0.011413521097787704; 0.027972308950302116; + 0.051727015600492421; 0.082502225484340941; 0.12007019910960293; + 0.16415283300752470; 0.21442376986779355; 0.27051082840644336; + 0.33199876341447887; 0.39843234186401943; 0.46931971407375483; + 0.54413605556657973; 0.62232745288031077; 0.70331500465597174; + 0.78649910768313447; 0.87126389619061517; 0.95698180152629142 |] - let a1 = a - 1.0 + // Weights for Gauss-Legendre quadrature + let gaussLegQuadWeights = [| + 0.0055657196642445571; 0.012915947284065419; 0.020181515297735382; + 0.027298621498568734; 0.034213810770299537; 0.040875750923643261; + 0.047235083490265582; 0.053244713977759692; 0.058860144245324798; + 0.064039797355015485; 0.068745323835736408; 0.072941885005653087; + 0.076598410645870640; 0.079687828912071670; 0.082187266704339706; + 0.084078218979661945; 0.085346685739338721; 0.085983275670394821 |] + + let ngau = gaussLegQuadWeights.Length + let a1 = a - 'T.One let lna1 = log a1 let sqrta1 = sqrt a1 - let gln = gammaLn a + let gln = _gammaLn a + let xu = - if (x > a1) then - max (a1 + 11.5*sqrta1) (x + 6.0*sqrta1) - else - max 0. ( min (a1 - 7.5*sqrta1) (x - 5.0*sqrta1)) + if x > a1 then + max (a1 + T<'T> 11.5 * sqrta1) (x + T<'T> 6.0 * sqrta1) + else + max 'T.Zero (min (a1 - T<'T> 7.5 * sqrta1) (x - T<'T> 5.0 * sqrta1)) - let rec gaussLegendre j sum = + let rec gaussLegendre j (sum: 'T) = if j < ngau then - let t = x + (xu - x) * gaussLegQuadY.[j] - let sum' = sum + gaussLegQuadWeights.[j] * exp(-(t - a1) + a1 * (log(t) - lna1)) - gaussLegendre (j+1) sum' + let t: 'T = x + (xu - x) * T<'T> gaussLegQuadY.[j] + let w: 'T = T<'T> gaussLegQuadWeights.[j] + let f: 'T = exp (-(t - a1) + a1 * (log t - lna1)) + gaussLegendre (j + 1) (sum + w * f) else - let ans = sum * (xu - x) * exp(a1 * (lna1 - 1.) - gln) - if psig then - if (ans > 0.0) then 1.0 - ans else -ans + let ans = sum * (xu - x) * exp (a1 * (lna1 - 'T.One) - gln) + if psig then + if ans > 'T.Zero then 'T.One - ans else -ans else - if (ans>=0.0) then ans else 1.0 + ans + if ans >= 'T.Zero then ans else 'T.One + ans - gaussLegendre 0 0.0 + gaussLegendre 0 'T.Zero /// Returns the incomplete gamma function P(a,X) (regularized gamma) // gammp -> GammaLowerIncomplete - let lowerIncompleteRegularized a x = - let ASWITCH=100. - //if (x < 0.0 || a <= 0.0) then failwith ("bad args in gammp") + let inline lowerIncompleteRegularized<'T when 'T :> Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : comparison > (a: 'T) (x: 'T) : 'T = + + let ASWITCH = T<'T> 100.0 match x with - | x when x = 0. -> 0. - | x when (x < 0.0 || a <= 0.0) -> nan - | x when (Ops.isPosInf x) -> 1. - | _ -> - if (x= 0.0) then 0.0 - elif (a >= ASWITCH) then + | _ when x = 'T.Zero -> 'T.Zero + | _ when x < 'T.Zero || a <= 'T.Zero -> T nan + | _ when 'T.IsPositiveInfinity x -> 'T.One + | _ -> + if a >= ASWITCH then gammpapprox a x true - - elif (x < a+1.0) then gser a x - else 1.0 - gcf a x + elif x < a + 'T.One then + gser a x + else + 'T.One - gcf a x /// Returns the incomplete gamma function Q(a,X) = 1 - P(a,X) (regularized gamma) // gammq -> GammaUpperIncomplete - let upperIncompleteRegularized a x = - let ASWITCH=100. - //if (x < 0.0 || a <= 0.0) then failwith ("bad args in gammp") + let inline upperIncompleteRegularized<'T when 'T :> Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : comparison > (a: 'T) ( x: 'T) : 'T = + + let ASWITCH = T<'T> 100.0 match x with - | x when x = 0. -> 1. - | x when (x < 0.0 || a <= 0.0) -> nan - | x when (Ops.isPosInf x) -> 0. - | _ -> - if (x= 0.0) then 1.0 - elif (a >= ASWITCH) then - 1.0 - gammpapprox a x true - - elif (x < a+1.0) then 1.0 - gser a x - else gcf a x + | _ when x = 'T.Zero -> 'T.One + | _ when x < 'T.Zero || a <= 'T.Zero -> T nan + | _ when 'T.IsPositiveInfinity x -> 'T.Zero + | _ -> + if a >= ASWITCH then + 'T.One - gammpapprox a x true + elif x < a + 'T.One then + 'T.One - gser a x + else + gcf a x /// Returns the lower incomplete gamma function
gamma(a,x) = int(exp(-t)t^(a-1),t=0..x) for real a > 0, x > 0. /// @@ -224,8 +324,14 @@ module Gamma = /// /// /// - let lowerIncomplete a x = - (lowerIncompleteRegularized a x) * gamma(a) + let inline lowerIncomplete<'T when 'T :> Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : comparison > (a: 'T) (x: 'T) : 'T = + + (lowerIncompleteRegularized a x) * (exp (_gammaLn a)) /// Returns the upper incomplete gamma function
Gamma(a,x) = int(exp(-t)t^(a-1),t=0..x) for real a > 0, x > 0. /// @@ -236,80 +342,121 @@ module Gamma = /// /// /// - let upperIncomplete a x = - (upperIncompleteRegularized a x) * gamma(a) + let inline upperIncomplete<'T when 'T :> Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : comparison > (a: 'T) (x: 'T) : 'T = + + (upperIncompleteRegularized a x) * (exp (_gammaLn a)) /// - /// Digamma function. + /// Optimized digamma function for x > 1. Skips reflection and tan-based correction. + /// Use when x is guaranteed to be greater than 1. /// - let rec digamma (x:float) = - match x with - | x when (x = 0.) -> - Double.NegativeInfinity - | x when (x < 0.) -> - digamma(1. - x) + Math.PI / Math.Tan(-Math.PI * x) - | _ -> - let nz = - - let q = x - let p = floor q + /// The input value (must be > 1). + /// The value of the digamma function ψ(x). + let inline digammaPositive<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T : comparison > (x: 'T) : 'T = + + let pi = 'T.Pi + let half = T<'T> 0.5 + let one = 'T.One + let zero = 'T.Zero + let ln = log<'T> + + let polyCoeffs = + [| + T 8.33333333333333333333E-2 + T -2.10927960927960927961E-2 + T 7.57575757575757575758E-3 + T -4.16666666666666666667E-3 + T 3.96825396825396825397E-3 + T -8.33333333333333333333E-3 + T 8.33333333333333333333E-2 + |] - if (p = q) then - raise (OverflowException("Function computation resulted in arithmetic overflow.")) + let ten = T<'T> 10.0 + let limit = T<'T> 1.0e17 - let nz = q - p - if (nz <> 0.5) then - if (nz > 0.5) then - let p = p + 1.0 - Math.PI / Math.Tan(System.Math.PI * (q - p)) - else - Math.PI / Math.Tan(System.Math.PI * nz) - else - 0.0 + let mutable s = x + let mutable w = zero - let x' = - if x <= 0 then - 1.0 - x - else - x + while s < ten do + w <- w + one / s + s <- s + one - let y = + if s < limit then + let z = one / (s * s) + let polv:'T = polyCoeffs |> Array.fold (fun acc c -> acc * z + c) zero + ln s - half / s - z * polv - w + else + ln s - half / s - w - if (x' <= 10.0 && x' = floor x') then - let n = floor x' - [1. .. (n-1.)] - |> Seq.fold (fun w y -> y + 1.0 / w) 0. - |> fun y -> y - 0.57721566490153286061 - else + /// + /// Full digamma function ψ(x), + /// + let inline digamma<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and Numerics.ITrigonometricFunctions<'T> + and 'T : comparison > (x: 'T) : 'T = + + let one = 'T.One + let zero = 'T.Zero + let half = T<'T> 0.5 + let ten = T<'T> 10.0 + let limit = T<'T> 1e17 + let pi = 'T.Pi + + let polyCoeffs = + [| + T 8.33333333333333333333E-2 + T -2.10927960927960927961E-2 + T 7.57575757575757575758E-3 + T -4.16666666666666666667E-3 + T 3.96825396825396825397E-3 + T -8.33333333333333333333E-3 + T 8.33333333333333333333E-2 + |] - let rec loop s w = - if (s < 10.0) then - let w = w + 1.0 / s - let s = s + 1.0 - loop s w - - else - if (s < 1.0E17) then - let z = 1.0 / (s * s) - let polv = - ((((((8.33333333333333333333E-2 - * z - 2.10927960927960927961E-2) - * z + 7.57575757575757575758E-3) - * z - 4.16666666666666666667E-3) - * z + 3.96825396825396825397E-3) - * z - 8.33333333333333333333E-3) - * z + 8.33333333333333333333E-2) - Math.Log(s) - 0.5 / s - (z * polv) - w - - else - Math.Log(s) - 0.5 / s - w - - loop x' 0. - - if (x < 0.) then - y - nz + let input = + if x < zero then + one - x // reflect else - y + x + + // Accumulate series if input is small + let mutable s = input + let mutable w = zero + while s < ten do + w <- w + one / s + s <- s + one + + // Asymptotic expansion + let result = + if s < limit then + let z = one / (s * s) + let polv: 'T = polyCoeffs |> Array.fold (fun acc c -> acc * z + c) zero + log<'T> s - half / s - z * polv - w + else + log<'T> s - half / s - w + + if x = zero then + T System.Double.NegativeInfinity + elif x < zero then + result + pi / 'T.Tan (-pi * x) + else + result /// @@ -322,7 +469,7 @@ module Gamma = /// made public under the GNU LGPL license. /// let rec trigamma (x:float) = - + //TODO: adapt to 'T match x with | x when (x < 0.) -> let v = Math.PI / Math.Sin(-Math.PI * x) From 1c82ab7d894ed76c2c8b78ebd4738764848a7bcd Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Mon, 23 Jun 2025 16:50:47 +0200 Subject: [PATCH 02/12] Update Beta.fs Marked _betaLn, _beta, betaLn, and beta functions as 'inline' to enable generic numeric support. --- src/FSharp.Stats/SpecialFunctions/Beta.fs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/FSharp.Stats/SpecialFunctions/Beta.fs b/src/FSharp.Stats/SpecialFunctions/Beta.fs index ab83be60..42774a9b 100644 --- a/src/FSharp.Stats/SpecialFunctions/Beta.fs +++ b/src/FSharp.Stats/SpecialFunctions/Beta.fs @@ -19,7 +19,7 @@ module Beta = /// /// The function input for approximating ln(B(z, w)) /// The function input for approximating ln(B(z, w)) - let _betaLn z w = (Gamma._gammaLn z) + (Gamma._gammaLn w) - (Gamma._gammaLn (z+w)) + let inline _betaLn (z: 'T) (w: 'T) = (Gamma._gammaLn z) + (Gamma._gammaLn w) - (Gamma._gammaLn (z+w)) /// /// Computes an approximation of the real value of the beta function using approximations for the gamma function using Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -29,7 +29,7 @@ module Beta = /// /// The function input for approximating B(z, w) /// The function input for approximating B(z, w) - let _beta z w = exp (_betaLn z w) + let inline _beta (z: 'T) (w: 'T) = exp (_betaLn z w) /// /// Computes an approximation of the real value of the log beta function using approximations for the gamma function using Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -40,7 +40,7 @@ module Beta = /// /// The function input for approximating ln(B(z, w)) /// The function input for approximating ln(B(z, w)) - let betaLn z w = (Gamma.gammaLn z) + (Gamma.gammaLn w) - (Gamma.gammaLn (z+w)) + let inline betaLn (z: 'T) (w: 'T) = (Gamma.gammaLn z) + (Gamma.gammaLn w) - (Gamma.gammaLn (z+w)) /// /// Computes an approximation of the real value of the beta function using approximations for the gamma function using Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -51,7 +51,7 @@ module Beta = /// /// The function input for approximating B(z, w) /// The function input for approximating B(z, w) - let beta z w = exp (betaLn z w) + let inline beta (z: 'T) (w: 'T) = exp (betaLn z w) // incomplete beta function /// From dc849cce84c05bee0ff0c500c0b9b8c32b012489 Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Tue, 24 Jun 2025 16:18:50 +0200 Subject: [PATCH 03/12] Start Update for Beta.lowerIncompleteRegularized<'T> --- src/FSharp.Stats/SpecialFunctions/Beta.fs | 188 ++++++++++++++++------ 1 file changed, 139 insertions(+), 49 deletions(-) diff --git a/src/FSharp.Stats/SpecialFunctions/Beta.fs b/src/FSharp.Stats/SpecialFunctions/Beta.fs index 42774a9b..2195885b 100644 --- a/src/FSharp.Stats/SpecialFunctions/Beta.fs +++ b/src/FSharp.Stats/SpecialFunctions/Beta.fs @@ -2,14 +2,15 @@ open System open FSharp.Stats +open FSharp.Stats.GenericMath /// The beta function B(p,q), or the beta integral (also called the Eulerian integral of the first kind) is defined by /// /// B(p, q) = (Γ(p) * Γ(q)) / Γ(p+q) module Beta = - let private EPS = 3.0e-8 // Precision.DoublePrecision; - let private FPMIN = 1.0e-30 // 0.0.Increment()/eps + let private EPS : 'T = T 3.0e-8 // Precision.DoublePrecision; + let private FPMIN : 'T = T 1.0e-30 // 0.0.Increment()/eps /// /// Computes an approximation of the real value of the log beta function using approximations for the gamma function using Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -19,7 +20,19 @@ module Beta = /// /// The function input for approximating ln(B(z, w)) /// The function input for approximating ln(B(z, w)) - let inline _betaLn (z: 'T) (w: 'T) = (Gamma._gammaLn z) + (Gamma._gammaLn w) - (Gamma._gammaLn (z+w)) + let inline _betaLn<'T when 'T :> Numerics.INumber<'T> + and 'T : (new: unit -> 'T) + and 'T : struct + and 'T : equality + and 'T :> ValueType + and 'T :> System.Numerics.IFloatingPoint<'T> + and 'T :> System.Numerics.IExponentialFunctions<'T> + and 'T :> System.Numerics.ILogarithmicFunctions<'T> + and 'T :> System.Numerics.IRootFunctions<'T> + and 'T :> System.Numerics.IPowerFunctions<'T> + and 'T : comparison> + (z: 'T) (w: 'T) = + (Gamma._gammaLn z) + (Gamma._gammaLn w) - (Gamma._gammaLn (z+w)) /// /// Computes an approximation of the real value of the beta function using approximations for the gamma function using Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -29,7 +42,19 @@ module Beta = /// /// The function input for approximating B(z, w) /// The function input for approximating B(z, w) - let inline _beta (z: 'T) (w: 'T) = exp (_betaLn z w) + let inline _beta<'T when 'T :> Numerics.INumber<'T> + and 'T : (new: unit -> 'T) + and 'T : struct + and 'T : equality + and 'T :> ValueType + and 'T :> System.Numerics.IFloatingPoint<'T> + and 'T :> System.Numerics.IExponentialFunctions<'T> + and 'T :> System.Numerics.ILogarithmicFunctions<'T> + and 'T :> System.Numerics.IRootFunctions<'T> + and 'T :> System.Numerics.IPowerFunctions<'T> + and 'T : comparison> + (z: 'T) (w: 'T) = + exp (_betaLn z w) /// /// Computes an approximation of the real value of the log beta function using approximations for the gamma function using Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -40,7 +65,19 @@ module Beta = /// /// The function input for approximating ln(B(z, w)) /// The function input for approximating ln(B(z, w)) - let inline betaLn (z: 'T) (w: 'T) = (Gamma.gammaLn z) + (Gamma.gammaLn w) - (Gamma.gammaLn (z+w)) + let inline betaLn<'T when 'T :> Numerics.INumber<'T> + and 'T : (new: unit -> 'T) + and 'T : struct + and 'T : equality + and 'T :> ValueType + and 'T :> System.Numerics.IFloatingPoint<'T> + and 'T :> System.Numerics.IExponentialFunctions<'T> + and 'T :> System.Numerics.ILogarithmicFunctions<'T> + and 'T :> System.Numerics.IRootFunctions<'T> + and 'T :> System.Numerics.IPowerFunctions<'T> + and 'T : comparison> + (z: 'T) (w: 'T) = + (Gamma.gammaLn z) + (Gamma.gammaLn w) - (Gamma.gammaLn (z+w)) /// /// Computes an approximation of the real value of the beta function using approximations for the gamma function using Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -51,7 +88,19 @@ module Beta = /// /// The function input for approximating B(z, w) /// The function input for approximating B(z, w) - let inline beta (z: 'T) (w: 'T) = exp (betaLn z w) + let inline beta<'T when 'T :> Numerics.INumber<'T> + and 'T : (new: unit -> 'T) + and 'T : struct + and 'T : equality + and 'T :> ValueType + and 'T :> System.Numerics.IFloatingPoint<'T> + and 'T :> System.Numerics.IExponentialFunctions<'T> + and 'T :> System.Numerics.ILogarithmicFunctions<'T> + and 'T :> System.Numerics.IRootFunctions<'T> + and 'T :> System.Numerics.IPowerFunctions<'T> + and 'T : comparison> + (z: 'T) (w: 'T) = + exp (betaLn z w) // incomplete beta function /// @@ -60,59 +109,100 @@ module Beta = /// The first Beta parameter, a positive real number. /// The second Beta parameter, a positive real number. /// The upper limit of the integral. - let lowerIncompleteRegularized a b x = - if (a < 0.0) then invalidArg "a" "Argument must not be negative" - if (b < 0.0) then invalidArg "b" "Argument must not be negative" - if (x < 0.0 || x > 1.0) then invalidArg "x" "Argument XY interval is inclusive" - let bt = - if (x = 0.0 || x = 1.0) then - 0.0 + let inline lowerIncompleteRegularized + (a: ^T when ^T :> Numerics.INumber<^T> + and ^T : (new: unit -> ^T) + and ^T : struct + and ^T : equality + and ^T :> ValueType + and ^T :> System.Numerics.IFloatingPoint<^T> + and ^T :> System.Numerics.IExponentialFunctions<^T> + and ^T :> System.Numerics.ILogarithmicFunctions<^T> + and ^T :> System.Numerics.IRootFunctions<^T> + and ^T :> System.Numerics.IPowerFunctions<^T> + and ^T : comparison) + (b: ^T when ^T :> Numerics.INumber<^T> + and ^T : (new: unit -> ^T) + and ^T : struct + and ^T : equality + and ^T :> ValueType + and ^T :> System.Numerics.IFloatingPoint<^T> + and ^T :> System.Numerics.IExponentialFunctions<^T> + and ^T :> System.Numerics.ILogarithmicFunctions<^T> + and ^T :> System.Numerics.IRootFunctions<^T> + and ^T :> System.Numerics.IPowerFunctions<^T> + and ^T : comparison) + (x: ^T when ^T :> Numerics.INumber<^T> + and ^T : (new: unit -> ^T) + and ^T : struct + and ^T : equality + and ^T :> ValueType + and ^T :> System.Numerics.IFloatingPoint<^T> + and ^T :> System.Numerics.IExponentialFunctions<^T> + and ^T :> System.Numerics.ILogarithmicFunctions<^T> + and ^T :> System.Numerics.IRootFunctions<^T> + and ^T :> System.Numerics.IPowerFunctions<^T> + and ^T : comparison) = + if (a < T 0.0) then invalidArg "a" "Argument must not be negative" + if (b < T 0.0) then invalidArg "b" "Argument must not be negative" + if (x < T 0.0 || x > T 1.0) then invalidArg "x" "Argument XY interval is inclusive" + let bt: 'T = + if (x = T 0.0 || x = T 1.0) then + T 0.0 else - exp (Gamma._gammaLn (a + b) - Gamma._gammaLn a - Gamma._gammaLn b + (a*Math.Log(x)) + (b*Math.Log(1.0 - x))) + exp (Gamma._gammaLn (a + b) - Gamma._gammaLn a - Gamma._gammaLn b + (a*log(x)) + (b*log(T 1.0 - x))) - let isSymmetryTransformation = ( x >= (a + 1.0)/(a + b + 2.0)) + let isSymmetryTransformation = ( x >= (a + T 1.0)/(a + b + T 2.0)) - let symmetryTransformation a b x = + let inline symmetryTransformation + (a: 'T) (b: 'T) (x: 'T) = let qab = a + b - let qap = a + 1.0 - let qam = a - 1.0 - let c = 1.0 - let d = - let tmp = 1.0 - (qab * x / qap) - if (abs tmp < FPMIN) then 1. / FPMIN else 1. / tmp + let qap = a + T 1.0 + let qam = a - T 1.0 + let c = T 1.0 + let d: 'T = + let tmp = T 1.0 - (qab * x / qap) + if (abs tmp < FPMIN) then T 1. / FPMIN else T 1. / tmp let h = d - let rec loop m mm d h c = - let aa = float m * (b - float m)*x/((qam + mm)*(a + mm)) + let rec loop (m: 'T) (mm: 'T) (d: 'T) (h: 'T) (c: 'T) = + let aa = m * (b - m)*x/((qam + mm)*(a + mm)) let d' = - let tmp = 1.0 + (aa*d) - if (abs tmp < FPMIN) then 1. / FPMIN else 1. / tmp + let tmp = T 1.0 + (aa*d) + if (abs tmp < FPMIN) then T 1. / FPMIN else T 1. / tmp let c' = - let tmp = 1.0 + (aa/c) + let tmp = T 1.0 + (aa/c) if (abs tmp < FPMIN) then FPMIN else tmp - let h' = h * d' * c' - let aa' = -(a + float m)*(qab + float m)*x/((a + mm)*(qap + mm)) - let d'' = - let tmp = 1.0 + (aa' * d') - if (abs tmp < FPMIN) then 1. / FPMIN else 1. / tmp - let c'' = - let tmp = 1.0 + (aa'/c') + let h': 'T = h * d' * c' + let aa': 'T = -(a + m)*(qab + m)*x/((a + mm)*(qap + mm)) + let d'': 'T = + let tmp = T 1.0 + (aa' * d') + if (abs tmp < FPMIN) then T 1. / FPMIN else T 1. / tmp + let c'': 'T = + let tmp = T 1.0 + (aa'/c') if (abs tmp < FPMIN) then FPMIN else tmp - let del = d''*c'' - let h'' = h' * del + let del: 'T = d''*c'' + let h'': 'T = h' * del - if abs (del - 1.0) <= EPS then - if isSymmetryTransformation then 1.0 - (bt*h''/a) else bt*h''/a + let tepmI: 'T = + bt*h''/a + + if abs (del - T 1.0) <= EPS then + + if isSymmetryTransformation then + T 1.0 - tepmI + else + tepmI else - if m < 140 then - loop (m+1) (mm+2.) d'' h'' c'' + if m < T 140 then + loop (m+T 1) (mm+T 2.) d'' h'' c'' else - if isSymmetryTransformation then 1.0 - (bt*h''/a) else bt*h''/a + if isSymmetryTransformation then T 1.0 - tepmI else tepmI - loop 1 2. d h c + loop (T 1) (T 2.) d h c if isSymmetryTransformation then - symmetryTransformation b a (1.0-x) + symmetryTransformation b a (T 1.0-x) else symmetryTransformation a b x @@ -123,7 +213,7 @@ module Beta = /// The first Beta parameter, a positive real number. /// The second Beta parameter, a positive real number. /// The upper limit of the integral. - let lowerIncomplete(a, b, x) = + let inline lowerIncomplete(a, b, x) : 'T= (lowerIncompleteRegularized a b x) * (beta a b) @@ -131,11 +221,11 @@ module Beta = /// Power series for incomplete beta integral. Use when b*x /// is small and x not too close to 1. /// - let powerSeries a b x = - let ai = 1.0 / a - let ui = (1.0 - b) * x - let t1 = ui / (a + 1.0) - let z = Ops.epsilon * ai + let inline powerSeries a b x = + let ai = T 1.0 / a + let ui = (T 1.0 - b) * x + let t1 = ui / (a + T 1.0) + let z = Ops.epsilon * ai |> T let rec loop u t v s n = if (abs v > z) then From bc3f0b7423b5b701c8c54136987911d284973c5b Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Tue, 24 Jun 2025 16:19:03 +0200 Subject: [PATCH 04/12] Start Update for TestStatistics for 'T --- .../Distributions/Continuous/StudentT.fs | 9 +++-- src/FSharp.Stats/Testing/TestStatistics.fs | 40 +++++++++++-------- 2 files changed, 29 insertions(+), 20 deletions(-) diff --git a/src/FSharp.Stats/Distributions/Continuous/StudentT.fs b/src/FSharp.Stats/Distributions/Continuous/StudentT.fs index d3053739..f6fa8b93 100644 --- a/src/FSharp.Stats/Distributions/Continuous/StudentT.fs +++ b/src/FSharp.Stats/Distributions/Continuous/StudentT.fs @@ -4,6 +4,7 @@ open System open FSharp.Stats open FSharp.Stats.Distributions open FSharp.Stats.Ops +open FSharp.Stats.GenericMath // ###### // Student's T-distribution @@ -17,10 +18,10 @@ open FSharp.Stats.Ops type StudentT = // Student's T-distribution helper functions. - static member CheckParam mu tau dof = - if System.Double.IsNaN(mu) || tau < 0.0 || System.Double.IsNaN(dof) || dof < 0. then + static member CheckParam (mu: 'T) (tau: 'T) (dof: 'T) = + if isNan(mu) || tau < T 0.0 || isNan(dof) || dof < T 0. then failwith "Student's T-distribution should be parametrized by mu, tau and dof > 0.0." - + /// Computes the mode. /// /// @@ -63,7 +64,7 @@ type StudentT = StudentT.CheckParam mu tau dof match dof with | df when System.Double.IsPositiveInfinity(df) -> tau*tau - | df when df > 2.0 -> dof*tau*tau/(dof-2.0) + | df when df > T 2.0 -> dof*tau*tau/(dof-T 2.0) | _ -> System.Double.PositiveInfinity /// Computes the standard deviation. diff --git a/src/FSharp.Stats/Testing/TestStatistics.fs b/src/FSharp.Stats/Testing/TestStatistics.fs index 955c125f..b6e4fddd 100644 --- a/src/FSharp.Stats/Testing/TestStatistics.fs +++ b/src/FSharp.Stats/Testing/TestStatistics.fs @@ -1,4 +1,6 @@ namespace FSharp.Stats.Testing +open System +open FSharp.Stats.GenericMath module TestStatistics = @@ -9,25 +11,31 @@ module TestStatistics = /// Creates a new T-Test for a given statistic /// with given degrees of freedom. /// - type TTestStatistics = { - /// The test statistic. - Statistic : float - /// The degrees of freedom for the numerator. - DegreesOfFreedom : float - /// One Tailed/Sided. - PValueLeft : float - /// One Tailed/Sided. - PValueRight : float - /// Two Tailed/Sided. - PValue : float - } + /// + /// The test statistic. + /// The degrees of freedom for the numerator. + /// One Tailed/Sided. + /// One Tailed/Sided. + /// Two Tailed/Sided. + type TTestStatistics<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T>> = + { + Statistic : 'T + DegreesOfFreedom : 'T + PValueLeft : 'T + PValueRight : 'T + PValue : 'T + } - let createTTest statistic dof = - let cdf = Distributions.Continuous.StudentT.CDF 0. 1. dof statistic - let pvalue = if statistic > 0. then 1. - cdf else cdf + let createTTest (statistic: 'T) (dof: 'T) = + let cdf: 'T = Distributions.Continuous.StudentT.CDF 0. 1. ('T.One dof) statistic + let pvalue: 'T = if statistic > T 0. then T 1. - cdf else cdf + {Statistic=statistic; DegreesOfFreedom=dof; PValueLeft=1. - pvalue; PValueRight=pvalue; PValue=pvalue*2.;} - /// /// Creates a new F-Test for a given statistic /// with given degrees of freedom. From eb3b9d281637a47b5cdac5d9fa01e0105286d318 Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Thu, 26 Jun 2025 16:14:51 +0200 Subject: [PATCH 05/12] Refactor to support generic numeric type 'T --- src/FSharp.Stats/SpecialFunctions/Beta.fs | 140 ++++++++-------------- 1 file changed, 52 insertions(+), 88 deletions(-) diff --git a/src/FSharp.Stats/SpecialFunctions/Beta.fs b/src/FSharp.Stats/SpecialFunctions/Beta.fs index 2195885b..480a879a 100644 --- a/src/FSharp.Stats/SpecialFunctions/Beta.fs +++ b/src/FSharp.Stats/SpecialFunctions/Beta.fs @@ -4,13 +4,18 @@ open System open FSharp.Stats open FSharp.Stats.GenericMath +// TODO: Generalize lowerIncompleteRegularized, lowerIncomplete, and powerSeries to support generic type 'T instead of using float + /// The beta function B(p,q), or the beta integral (also called the Eulerian integral of the first kind) is defined by /// /// B(p, q) = (Γ(p) * Γ(q)) / Γ(p+q) module Beta = - let private EPS : 'T = T 3.0e-8 // Precision.DoublePrecision; - let private FPMIN : 'T = T 1.0e-30 // 0.0.Increment()/eps + // let private EPS : 'T = T 3.0e-8 // Precision.DoublePrecision; + // let private FPMIN : 'T = T 1.0e-30 // 0.0.Increment()/eps + + let private f_EPS = T 3.0e-8 // Precision.DoublePrecision; + let private f_FPMIN = T 1.0e-30 // 0.0.Increment()/eps /// /// Computes an approximation of the real value of the log beta function using approximations for the gamma function using Lanczos Coefficients described in Numerical Recipes (Press et al) @@ -109,100 +114,59 @@ module Beta = /// The first Beta parameter, a positive real number. /// The second Beta parameter, a positive real number. /// The upper limit of the integral. - let inline lowerIncompleteRegularized - (a: ^T when ^T :> Numerics.INumber<^T> - and ^T : (new: unit -> ^T) - and ^T : struct - and ^T : equality - and ^T :> ValueType - and ^T :> System.Numerics.IFloatingPoint<^T> - and ^T :> System.Numerics.IExponentialFunctions<^T> - and ^T :> System.Numerics.ILogarithmicFunctions<^T> - and ^T :> System.Numerics.IRootFunctions<^T> - and ^T :> System.Numerics.IPowerFunctions<^T> - and ^T : comparison) - (b: ^T when ^T :> Numerics.INumber<^T> - and ^T : (new: unit -> ^T) - and ^T : struct - and ^T : equality - and ^T :> ValueType - and ^T :> System.Numerics.IFloatingPoint<^T> - and ^T :> System.Numerics.IExponentialFunctions<^T> - and ^T :> System.Numerics.ILogarithmicFunctions<^T> - and ^T :> System.Numerics.IRootFunctions<^T> - and ^T :> System.Numerics.IPowerFunctions<^T> - and ^T : comparison) - (x: ^T when ^T :> Numerics.INumber<^T> - and ^T : (new: unit -> ^T) - and ^T : struct - and ^T : equality - and ^T :> ValueType - and ^T :> System.Numerics.IFloatingPoint<^T> - and ^T :> System.Numerics.IExponentialFunctions<^T> - and ^T :> System.Numerics.ILogarithmicFunctions<^T> - and ^T :> System.Numerics.IRootFunctions<^T> - and ^T :> System.Numerics.IPowerFunctions<^T> - and ^T : comparison) = - if (a < T 0.0) then invalidArg "a" "Argument must not be negative" - if (b < T 0.0) then invalidArg "b" "Argument must not be negative" - if (x < T 0.0 || x > T 1.0) then invalidArg "x" "Argument XY interval is inclusive" - let bt: 'T = - if (x = T 0.0 || x = T 1.0) then - T 0.0 + let lowerIncompleteRegularized a b x = + if (a < 0.0) then invalidArg "a" "Argument must not be negative" + if (b < 0.0) then invalidArg "b" "Argument must not be negative" + if (x < 0.0 || x > 1.0) then invalidArg "x" "Argument XY interval is inclusive" + let bt = + if (x = 0.0 || x = 1.0) then + 0.0 else - exp (Gamma._gammaLn (a + b) - Gamma._gammaLn a - Gamma._gammaLn b + (a*log(x)) + (b*log(T 1.0 - x))) + exp (Gamma._gammaLn (a + b) - Gamma._gammaLn a - Gamma._gammaLn b + (a*Math.Log(x)) + (b*Math.Log(1.0 - x))) - let isSymmetryTransformation = ( x >= (a + T 1.0)/(a + b + T 2.0)) + let isSymmetryTransformation = ( x >= (a + 1.0)/(a + b + 2.0)) - let inline symmetryTransformation - (a: 'T) (b: 'T) (x: 'T) = + let symmetryTransformation a b x = let qab = a + b - let qap = a + T 1.0 - let qam = a - T 1.0 - let c = T 1.0 - let d: 'T = - let tmp = T 1.0 - (qab * x / qap) - if (abs tmp < FPMIN) then T 1. / FPMIN else T 1. / tmp + let qap = a + 1.0 + let qam = a - 1.0 + let c = 1.0 + let d = + let tmp = 1.0 - (qab * x / qap) + if (abs tmp < f_FPMIN) then 1. / f_FPMIN else 1. / tmp let h = d - let rec loop (m: 'T) (mm: 'T) (d: 'T) (h: 'T) (c: 'T) = - let aa = m * (b - m)*x/((qam + mm)*(a + mm)) + let rec loop m mm d h c = + let aa = float m * (b - float m)*x/((qam + mm)*(a + mm)) let d' = - let tmp = T 1.0 + (aa*d) - if (abs tmp < FPMIN) then T 1. / FPMIN else T 1. / tmp + let tmp = 1.0 + (aa*d) + if (abs tmp < f_FPMIN) then 1. / f_FPMIN else 1. / tmp let c' = - let tmp = T 1.0 + (aa/c) - if (abs tmp < FPMIN) then FPMIN else tmp - let h': 'T = h * d' * c' - let aa': 'T = -(a + m)*(qab + m)*x/((a + mm)*(qap + mm)) - let d'': 'T = - let tmp = T 1.0 + (aa' * d') - if (abs tmp < FPMIN) then T 1. / FPMIN else T 1. / tmp - let c'': 'T = - let tmp = T 1.0 + (aa'/c') - if (abs tmp < FPMIN) then FPMIN else tmp + let tmp = 1.0 + (aa/c) + if (abs tmp < f_FPMIN) then f_FPMIN else tmp + let h' = h * d' * c' + let aa' = -(a + float m)*(qab + float m)*x/((a + mm)*(qap + mm)) + let d'' = + let tmp = 1.0 + (aa' * d') + if (abs tmp < f_FPMIN) then 1. / f_FPMIN else 1. / tmp + let c'' = + let tmp = 1.0 + (aa'/c') + if (abs tmp < f_FPMIN) then f_FPMIN else tmp - let del: 'T = d''*c'' - let h'': 'T = h' * del + let del = d''*c'' + let h'' = h' * del - let tepmI: 'T = - bt*h''/a - - if abs (del - T 1.0) <= EPS then - - if isSymmetryTransformation then - T 1.0 - tepmI - else - tepmI + if abs (del - 1.0) <= f_EPS then + if isSymmetryTransformation then 1.0 - (bt*h''/a) else bt*h''/a else - if m < T 140 then - loop (m+T 1) (mm+T 2.) d'' h'' c'' + if m < 140 then + loop (m+1) (mm+2.) d'' h'' c'' else - if isSymmetryTransformation then T 1.0 - tepmI else tepmI + if isSymmetryTransformation then 1.0 - (bt*h''/a) else bt*h''/a - loop (T 1) (T 2.) d h c + loop 1 2. d h c if isSymmetryTransformation then - symmetryTransformation b a (T 1.0-x) + symmetryTransformation b a (1.0-x) else symmetryTransformation a b x @@ -213,7 +177,7 @@ module Beta = /// The first Beta parameter, a positive real number. /// The second Beta parameter, a positive real number. /// The upper limit of the integral. - let inline lowerIncomplete(a, b, x) : 'T= + let lowerIncomplete(a, b, x) = (lowerIncompleteRegularized a b x) * (beta a b) @@ -221,11 +185,11 @@ module Beta = /// Power series for incomplete beta integral. Use when b*x /// is small and x not too close to 1. /// - let inline powerSeries a b x = - let ai = T 1.0 / a - let ui = (T 1.0 - b) * x - let t1 = ui / (a + T 1.0) - let z = Ops.epsilon * ai |> T + let powerSeries a b x = + let ai = 1.0 / a + let ui = (1.0 - b) * x + let t1 = ui / (a + 1.0) + let z = Ops.epsilon * ai let rec loop u t v s n = if (abs v > z) then From 64ac64ae6fda326d23f04819ce2e596181fc3a1f Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Thu, 26 Jun 2025 16:15:31 +0200 Subject: [PATCH 06/12] Partial update to support generic numeric type 'T in TestStatistics --- .../Distributions/Continuous/StudentT.fs | 7 +- src/FSharp.Stats/Testing/TestStatistics.fs | 184 +++++++++++++----- 2 files changed, 142 insertions(+), 49 deletions(-) diff --git a/src/FSharp.Stats/Distributions/Continuous/StudentT.fs b/src/FSharp.Stats/Distributions/Continuous/StudentT.fs index f6fa8b93..dd5df8fb 100644 --- a/src/FSharp.Stats/Distributions/Continuous/StudentT.fs +++ b/src/FSharp.Stats/Distributions/Continuous/StudentT.fs @@ -4,7 +4,6 @@ open System open FSharp.Stats open FSharp.Stats.Distributions open FSharp.Stats.Ops -open FSharp.Stats.GenericMath // ###### // Student's T-distribution @@ -18,8 +17,8 @@ open FSharp.Stats.GenericMath type StudentT = // Student's T-distribution helper functions. - static member CheckParam (mu: 'T) (tau: 'T) (dof: 'T) = - if isNan(mu) || tau < T 0.0 || isNan(dof) || dof < T 0. then + static member CheckParam mu tau dof = + if isNan(mu) || tau < 0.0 || isNan(dof) || dof < 0. then failwith "Student's T-distribution should be parametrized by mu, tau and dof > 0.0." /// Computes the mode. @@ -64,7 +63,7 @@ type StudentT = StudentT.CheckParam mu tau dof match dof with | df when System.Double.IsPositiveInfinity(df) -> tau*tau - | df when df > T 2.0 -> dof*tau*tau/(dof-T 2.0) + | df when df > 2.0 -> dof*tau*tau/(dof-2.0) | _ -> System.Double.PositiveInfinity /// Computes the standard deviation. diff --git a/src/FSharp.Stats/Testing/TestStatistics.fs b/src/FSharp.Stats/Testing/TestStatistics.fs index b6e4fddd..d620c286 100644 --- a/src/FSharp.Stats/Testing/TestStatistics.fs +++ b/src/FSharp.Stats/Testing/TestStatistics.fs @@ -2,6 +2,7 @@ namespace FSharp.Stats.Testing open System open FSharp.Stats.GenericMath +// TODO: Update specific distributions to support generic type 'T to avoid explicit float casting module TestStatistics = @@ -21,7 +22,9 @@ module TestStatistics = and Numerics.IFloatingPoint<'T> and Numerics.IExponentialFunctions<'T> and Numerics.IRootFunctions<'T> - and Numerics.IPowerFunctions<'T>> = + and Numerics.IPowerFunctions<'T> + and 'T: (static member op_Explicit: ^T -> float) + and 'T : comparison > = { Statistic : 'T DegreesOfFreedom : 'T @@ -30,71 +33,162 @@ module TestStatistics = PValue : 'T } - let createTTest (statistic: 'T) (dof: 'T) = - let cdf: 'T = Distributions.Continuous.StudentT.CDF 0. 1. ('T.One dof) statistic - let pvalue: 'T = if statistic > T 0. then T 1. - cdf else cdf + let inline createTTest<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.ILogarithmicFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T: (static member op_Explicit: ^T -> float) + and 'T : comparison > + (statistic: 'T) (dof: 'T) = + let fstatistic = toFloat statistic + let fdof = toFloat dof + + let cdf = Distributions.Continuous.StudentT.CDF 0. 1. fdof fstatistic |> T + let pvalue = if fstatistic > 0. then 1. - cdf else cdf - {Statistic=statistic; DegreesOfFreedom=dof; PValueLeft=1. - pvalue; PValueRight=pvalue; PValue=pvalue*2.;} + { + Statistic=statistic; + DegreesOfFreedom=dof; + PValueLeft= T(1. - pvalue); + PValueRight=T(pvalue); + PValue=T(pvalue*2.); + } /// /// Creates a new F-Test for a given statistic /// with given degrees of freedom. /// - type FTestStatistics = { - /// The test statistic. - Statistic : float - /// The degrees of freedom for the numerator. - DegreesOfFreedom1 : float - /// The degrees of freedom for the denominator. - DegreesOfFreedom2 : float - PValue : float - PValueTwoTailed : float - } + /// + /// The test statistic. + /// The degrees of freedom for the numerator. + /// The degrees of freedom for the denominator. + type FTestStatistics<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T: (static member op_Explicit: ^T -> float) + and 'T : comparison > = + { + Statistic : 'T + DegreesOfFreedom1 : 'T + DegreesOfFreedom2 : 'T + PValue : 'T + PValueTwoTailed : 'T + } - let createFTest statistic dof1 dof2 = - let cdf = Distributions.Continuous.F.CDF dof1 dof2 statistic + let inline createFTest<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T: (static member op_Explicit: ^T -> float) + and 'T : comparison > + (statistic: 'T) (dof1: 'T) (dof2: 'T) = + let fstatistic = toFloat statistic + let fdof1 = toFloat dof1 + let fdof2 = toFloat dof2 + let cdf = Distributions.Continuous.F.CDF fdof1 fdof2 fstatistic let pvalue = 1. - cdf let pvalueTwoTailed = pvalue * 2. - {Statistic=statistic; DegreesOfFreedom1=dof1; DegreesOfFreedom2=dof2; PValue=pvalue; PValueTwoTailed = pvalueTwoTailed} + { + Statistic=statistic; + DegreesOfFreedom1=dof1; + DegreesOfFreedom2=dof2; + PValue=T pvalue; + PValueTwoTailed = T pvalueTwoTailed + } /// /// Computes the Chi-Square test statistics for a given statistic /// with given degrees of freedom. /// - type ChiSquareStatistics = { - /// The test statistic. - Statistic : float - /// The degrees of freedom for the numerator. - DegreesOfFreedom : float - /// One Tailed/Sided. - PValueLeft : float - /// One Tailed/Sided. - PValueRight : float - /// Two Tailed/Sided. - PValue : float - } + /// + /// The test statistic. + /// The degrees of freedom for the numerator. + /// One Tailed/Sided. + /// One Tailed/Sided. + /// Two Tailed/Sided. + type ChiSquareStatistics<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T: (static member op_Explicit: ^T -> float) + and 'T : comparison > = + { + Statistic : 'T + DegreesOfFreedom : 'T + /// one tailed/sided chiSquare pValue + PValueLeft : 'T + /// one tailed/sided chiSquare pValue (default) + PValueRight : 'T + /// two tailed/sided chiSquare pValue + PValue : 'T + } - let createChiSquare statistic dof = - let cdf = Distributions.Continuous.ChiSquared.CDF dof statistic - let pvalue = if statistic > 0. then 1. - cdf else cdf - {Statistic = statistic; DegreesOfFreedom = dof; PValueLeft = 1. - pvalue; PValueRight = pvalue; PValue = pvalue * 2.} + let inline createChiSquare<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T: (static member op_Explicit: ^T -> float) + and 'T : comparison > + (statistic:'T) (dof:'T) = + let fstatistic = toFloat statistic + let fdof = toFloat dof + + + let cdf = Distributions.Continuous.ChiSquared.CDF fdof fstatistic + let pvalue = if fstatistic > 0. then 1. - cdf else cdf + { + Statistic = statistic; + DegreesOfFreedom = dof; + PValueLeft =T(1. - pvalue); + PValueRight = T pvalue; + PValue = T(pvalue * 2.) + } /// /// Computes the Wilcoxon test statistics for a given statistic. /// - type WilcoxonTestStatistics = { - /// The test statistic. - Statistic : float - PValueLeft : float - PValueRight : float - /// Two Tailed/Sided. - PValueTwoTailed : float - } - let createWilcoxon statistic = - let cdf = Distributions.Continuous.Normal.CDF 0. 1. statistic + /// The test statistic. + /// One Tailed/Sided. + /// Two Tailed/Sided. + type WilcoxonTestStatistics<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T: (static member op_Explicit: ^T -> float) + and 'T : comparison > = + { + Statistic : 'T + PValueLeft : 'T + PValueRight : 'T + PValueTwoTailed : 'T + } + let inline createWilcoxon<'T when 'T :> Numerics.INumber<'T> + and Numerics.IFloatingPoint<'T> + and Numerics.IExponentialFunctions<'T> + and Numerics.IRootFunctions<'T> + and Numerics.IPowerFunctions<'T> + and 'T: (static member op_Explicit: ^T -> float) + and 'T : comparison > + (statistic:'T) = + let fstatistic = toFloat statistic + + let cdf = Distributions.Continuous.Normal.CDF 0. 1. fstatistic let pvalue = 1.- cdf let pvalueTwoTailed = pvalue * 2. - {Statistic=statistic; PValueLeft=pvalue;PValueRight = cdf; PValueTwoTailed = pvalueTwoTailed} \ No newline at end of file + { + Statistic=statistic; + PValueLeft=T pvalue; + PValueRight = T cdf; + PValueTwoTailed = T pvalueTwoTailed + } \ No newline at end of file From 8080cea2ec33b7eec5b3770617599cbb43448771 Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Thu, 26 Jun 2025 16:15:43 +0200 Subject: [PATCH 07/12] Update SpecialFunctions.fs --- tests/FSharp.Stats.Tests/SpecialFunctions.fs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/FSharp.Stats.Tests/SpecialFunctions.fs b/tests/FSharp.Stats.Tests/SpecialFunctions.fs index 1e1bd132..9f64f2d3 100644 --- a/tests/FSharp.Stats.Tests/SpecialFunctions.fs +++ b/tests/FSharp.Stats.Tests/SpecialFunctions.fs @@ -292,7 +292,7 @@ let factorialTests = Expect.floatClose Accuracy.high (Factorial._factorialLn 6942) 54467.727976695301612523565124699078303834231913072759124392135342 "factorialLn of large number failed" ) testCase "_ln(0!) = 0" (fun _ -> - Expect.equal (Factorial._factorialLn 0) 0. "Expected factorialLn of 0 to be 1." + Expect.floatClose Accuracy.high (Factorial._factorialLn 0) 0. "Expected factorialLn of 0 to be 1." ) testCase "_ln(69!)" (fun _ -> Expect.floatClose Accuracy.high 226.19054832372759333227016852232261788323276357495863628461257077 (Factorial._factorialLn 69) "Expected factorialLn of 69 to be 226.19054832372759333227016852232261788323276357495863628461257077" @@ -310,7 +310,7 @@ let factorialTests = Expect.floatClose Accuracy.high (Factorial.factorialLn 6942) 54467.727976695301612523565124699078303834231913072759124392135342 "factorialLn of large number failed" ) testCase "ln(0!) = 0" (fun _ -> - Expect.equal (Factorial.factorialLn 0) 0. "Expected factorialLn of 0 to be 1." + Expect.floatClose Accuracy.high (Factorial.factorialLn 0) 0. "Expected factorialLn of 0 to be 1." ) testCase "ln(69!)" (fun _ -> Expect.floatClose Accuracy.high 226.19054832372759333227016852232261788323276357495863628461257077 (Factorial.factorialLn 69) "Expected factorialLn of 69 to be 226.19054832372759333227016852232261788323276357495863628461257077" From 77ba43eec8d3c72dcda047b3c3e298a51be7f4af Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Mon, 20 Oct 2025 11:18:40 +0200 Subject: [PATCH 08/12] Update FSharp.Stats.GenericMath to FsMath.GenericMath --- src/FSharp.Stats/SpecialFunctions/Beta.fs | 3 ++- src/FSharp.Stats/SpecialFunctions/Gamma.fs | 3 ++- src/FSharp.Stats/Testing/TestStatistics.fs | 3 ++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/FSharp.Stats/SpecialFunctions/Beta.fs b/src/FSharp.Stats/SpecialFunctions/Beta.fs index 480a879a..4f88f0f2 100644 --- a/src/FSharp.Stats/SpecialFunctions/Beta.fs +++ b/src/FSharp.Stats/SpecialFunctions/Beta.fs @@ -2,7 +2,8 @@ open System open FSharp.Stats -open FSharp.Stats.GenericMath +open FsMath +open FsMath.GenericMath // TODO: Generalize lowerIncompleteRegularized, lowerIncomplete, and powerSeries to support generic type 'T instead of using float diff --git a/src/FSharp.Stats/SpecialFunctions/Gamma.fs b/src/FSharp.Stats/SpecialFunctions/Gamma.fs index 0e662773..724c9738 100644 --- a/src/FSharp.Stats/SpecialFunctions/Gamma.fs +++ b/src/FSharp.Stats/SpecialFunctions/Gamma.fs @@ -1,7 +1,8 @@ namespace FSharp.Stats.SpecialFunctions open System -open FSharp.Stats.GenericMath +open FsMath +open FsMath.GenericMath /// Approximations for the gamma function and related functions. /// diff --git a/src/FSharp.Stats/Testing/TestStatistics.fs b/src/FSharp.Stats/Testing/TestStatistics.fs index d620c286..2670a9fe 100644 --- a/src/FSharp.Stats/Testing/TestStatistics.fs +++ b/src/FSharp.Stats/Testing/TestStatistics.fs @@ -1,6 +1,7 @@ namespace FSharp.Stats.Testing open System -open FSharp.Stats.GenericMath +open FsMath +open FsMath.GenericMath // TODO: Update specific distributions to support generic type 'T to avoid explicit float casting From b956775b17e4489fbdea57e2c29077714fec897b Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Tue, 21 Oct 2025 09:57:30 +0200 Subject: [PATCH 09/12] Update FsMath Reference --- src/FSharp.Stats/FSharp.Stats.fsproj | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FSharp.Stats/FSharp.Stats.fsproj b/src/FSharp.Stats/FSharp.Stats.fsproj index c7dc2ee1..a994bcd4 100644 --- a/src/FSharp.Stats/FSharp.Stats.fsproj +++ b/src/FSharp.Stats/FSharp.Stats.fsproj @@ -172,7 +172,7 @@ - + From 47d93e90ab347ef85472b676e397a9e5f8092aec Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Tue, 21 Oct 2025 10:06:50 +0200 Subject: [PATCH 10/12] Update FsMath Reference in Docs --- docs/BasicStats.fsx | 2 +- docs/Clustering.fsx | 2 +- docs/ComparisonMetrics.fsx | 4 ++-- docs/Correlation.fsx | 6 +++--- docs/Covariance.fsx | 2 +- docs/CrossValidation.fsx | 2 +- docs/Differentiation.fsx | 2 +- docs/Distributions.fsx | 2 +- docs/Fitting.fsx | 2 +- docs/GoodnessOfFit.fsx | 2 +- docs/GrowthCurve.fsx | 2 +- docs/Imputation.fsx | 4 ++-- docs/Integration.fsx | 2 +- docs/Interpolation.fsx | 2 +- docs/Intervals.fsx | 2 +- docs/Normalization.fsx | 2 +- docs/Optimization.fsx | 2 +- docs/Quantiles.fsx | 2 +- docs/Rank.fsx | 2 +- docs/Signal.fsx | 2 +- docs/Testing.fsx | 2 +- docs/index.fsx | 2 +- 22 files changed, 26 insertions(+), 26 deletions(-) diff --git a/docs/BasicStats.fsx b/docs/BasicStats.fsx index fe1a11be..2f66fcca 100644 --- a/docs/BasicStats.fsx +++ b/docs/BasicStats.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" (*** condition: ipynb ***) diff --git a/docs/Clustering.fsx b/docs/Clustering.fsx index d846b6c6..e63ff823 100644 --- a/docs/Clustering.fsx +++ b/docs/Clustering.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/ComparisonMetrics.fsx b/docs/ComparisonMetrics.fsx index b4a141ae..056e7970 100644 --- a/docs/ComparisonMetrics.fsx +++ b/docs/ComparisonMetrics.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" @@ -27,7 +27,7 @@ Plotly.NET.Defaults.DefaultDisplayOptions <- #if IPYNB #r "nuget: Plotly.NET, 4.0.0" #r "nuget: Plotly.NET.Interactive, 4.0.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #r "nuget: FSharp.Stats" #endif // IPYNB diff --git a/docs/Correlation.fsx b/docs/Correlation.fsx index 22442c85..34a52d44 100644 --- a/docs/Correlation.fsx +++ b/docs/Correlation.fsx @@ -14,11 +14,11 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" Plotly.NET.Defaults.DefaultDisplayOptions <- @@ -29,7 +29,7 @@ Plotly.NET.Defaults.DefaultDisplayOptions <- #r "nuget: Plotly.NET, 4.0.0" #r "nuget: Plotly.NET.Interactive, 4.0.0" #r "nuget: FSharp.Stats" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #endif // IPYNB (** diff --git a/docs/Covariance.fsx b/docs/Covariance.fsx index 7d14e30c..74702d85 100644 --- a/docs/Covariance.fsx +++ b/docs/Covariance.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/CrossValidation.fsx b/docs/CrossValidation.fsx index c97a5e2b..7202514b 100644 --- a/docs/CrossValidation.fsx +++ b/docs/CrossValidation.fsx @@ -15,7 +15,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Differentiation.fsx b/docs/Differentiation.fsx index 6484e9a4..2b98ff47 100644 --- a/docs/Differentiation.fsx +++ b/docs/Differentiation.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Distributions.fsx b/docs/Distributions.fsx index 3f769378..120e9edd 100644 --- a/docs/Distributions.fsx +++ b/docs/Distributions.fsx @@ -15,7 +15,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #r "nuget: Plotly.NET, 4.0.0" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" diff --git a/docs/Fitting.fsx b/docs/Fitting.fsx index fcf5c6bb..cee5dc91 100644 --- a/docs/Fitting.fsx +++ b/docs/Fitting.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/GoodnessOfFit.fsx b/docs/GoodnessOfFit.fsx index 2668cbdf..005715b0 100644 --- a/docs/GoodnessOfFit.fsx +++ b/docs/GoodnessOfFit.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/GrowthCurve.fsx b/docs/GrowthCurve.fsx index e5ad0513..d1ecb4f8 100644 --- a/docs/GrowthCurve.fsx +++ b/docs/GrowthCurve.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Imputation.fsx b/docs/Imputation.fsx index 528f065e..3d28a74f 100644 --- a/docs/Imputation.fsx +++ b/docs/Imputation.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" @@ -28,7 +28,7 @@ Plotly.NET.Defaults.DefaultDisplayOptions <- #r "nuget: Plotly.NET, 4.0.0" #r "nuget: Plotly.NET.Interactive, 4.0.0" #r "nuget: FSharp.Stats" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" open FsMath open Plotly.NET diff --git a/docs/Integration.fsx b/docs/Integration.fsx index 78523e31..d8667a4c 100644 --- a/docs/Integration.fsx +++ b/docs/Integration.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Interpolation.fsx b/docs/Interpolation.fsx index 54bedc99..53f35af5 100644 --- a/docs/Interpolation.fsx +++ b/docs/Interpolation.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Intervals.fsx b/docs/Intervals.fsx index 51c2b1ca..217336a1 100644 --- a/docs/Intervals.fsx +++ b/docs/Intervals.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Normalization.fsx b/docs/Normalization.fsx index 4a2229d1..4c68a45b 100644 --- a/docs/Normalization.fsx +++ b/docs/Normalization.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Optimization.fsx b/docs/Optimization.fsx index 4d6589a2..e0eeb3bf 100644 --- a/docs/Optimization.fsx +++ b/docs/Optimization.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Quantiles.fsx b/docs/Quantiles.fsx index a175be6d..511d3fd0 100644 --- a/docs/Quantiles.fsx +++ b/docs/Quantiles.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Rank.fsx b/docs/Rank.fsx index 63ee0619..fdb1e7ae 100644 --- a/docs/Rank.fsx +++ b/docs/Rank.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Signal.fsx b/docs/Signal.fsx index ce3e7f2d..9da9cff7 100644 --- a/docs/Signal.fsx +++ b/docs/Signal.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/Testing.fsx b/docs/Testing.fsx index 2e33df7d..ebd363e8 100644 --- a/docs/Testing.fsx +++ b/docs/Testing.fsx @@ -14,7 +14,7 @@ categoryindex: 0 #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" diff --git a/docs/index.fsx b/docs/index.fsx index 5b6146af..c7ab85cd 100644 --- a/docs/index.fsx +++ b/docs/index.fsx @@ -5,7 +5,7 @@ #r "nuget: FSharpAux, 2.0.0" #r "nuget: FSharpAux.IO, 2.0.0" #r "nuget: OptimizedPriorityQueue, 5.1.0" -#r "nuget: FsMath, 0.0.1" +#r "nuget: FsMath, 0.0.2" #I "../src/FSharp.Stats/bin/Release/.net8.0/" #r "FSharp.Stats.dll" #r "nuget: Plotly.NET, 4.0.0" From 00d3a4039dd9781eb9f0c575d79aac282c9b06bf Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Tue, 21 Oct 2025 10:16:00 +0200 Subject: [PATCH 11/12] Fix incorrect function calls in docs --- docs/ComparisonMetrics.fsx | 2 +- docs/CrossValidation.fsx | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/ComparisonMetrics.fsx b/docs/ComparisonMetrics.fsx index 056e7970..67be4665 100644 --- a/docs/ComparisonMetrics.fsx +++ b/docs/ComparisonMetrics.fsx @@ -154,7 +154,7 @@ let mlcm = [2; 0; 4] ] |> array2D - |> Matrix.Generic.ofArray2D + |> Matrix.ofArray2D ) ) (** diff --git a/docs/CrossValidation.fsx b/docs/CrossValidation.fsx index 7202514b..64911c18 100644 --- a/docs/CrossValidation.fsx +++ b/docs/CrossValidation.fsx @@ -357,7 +357,7 @@ let shuffleAndSplitPolynomial p iterations (xData: Vector) (yData: Vector let getFitFuncPol xTrain yTrain (xTest:Vector) = getFitFuncPolynomial xTrain yTrain xTest order - CrossValidation.shuffelAndSplit p iterations xDataMat yData getFitFuncPol error Seq.stDev + CrossValidation.shuffleAndSplit p iterations xDataMat yData getFitFuncPol error Seq.stDev //creates an output for 5 iterations where random 20 % of the data set are taken as testing data set let sasPolynomial order = shuffleAndSplitPolynomial 0.2 5 xV yV order @@ -368,7 +368,7 @@ let shuffleAndSplitSpline p iterations (xData: Vector) (yData: Vector) = getFitFuncSpline xDat yDat xDatTrain lambda - CrossValidation.shuffelAndSplit p iterations xDataMat yData getFitFuncSpl errorSpl Seq.stDev + CrossValidation.shuffleAndSplit p iterations xDataMat yData getFitFuncSpl errorSpl Seq.stDev //creates an output for 5 iterations where random 20 % of the data set are taken as testing data set let sasSpline lambda = shuffleAndSplitSpline 0.2 5 xV yV lambda From 050e07a149c19dd86d3cea997ec4091ad1a86972 Mon Sep 17 00:00:00 2001 From: Christopher Lux <59606965+LibraChris@users.noreply.github.com> Date: Tue, 21 Oct 2025 11:18:01 +0200 Subject: [PATCH 12/12] Fix XML Comments --- .../Distributions/Discrete/Hypergeometric.fs | 5 +- .../Discrete/NegativeBinomial.fs | 7 +-- src/FSharp.Stats/Fitting/CrossValidation.fs | 1 + src/FSharp.Stats/Fitting/LinearRegression.fs | 1 + .../Fitting/LogisticRegression.fs | 1 + src/FSharp.Stats/ML/Imputation.fs | 5 ++ .../ML/Unsupervised/ClusterNumber.fs | 3 + src/FSharp.Stats/ML/Unsupervised/KNN.fs | 2 + src/FSharp.Stats/Testing/TestStatistics.fs | 62 ++++++++++++------- 9 files changed, 58 insertions(+), 29 deletions(-) diff --git a/src/FSharp.Stats/Distributions/Discrete/Hypergeometric.fs b/src/FSharp.Stats/Distributions/Discrete/Hypergeometric.fs index c6431971..cfda78d8 100644 --- a/src/FSharp.Stats/Distributions/Discrete/Hypergeometric.fs +++ b/src/FSharp.Stats/Distributions/Discrete/Hypergeometric.fs @@ -250,9 +250,8 @@ type Hypergeometric = static member ToString N K n = sprintf "Hypergeometric(N = %i, K = %i, n = %i)" N K n - /// Initializes a hypergeometric distribution. - /// - /// + /// + /// Initializes a hypergeometric distribution. /// The hypergeometric distribution is a discrete probability distribution /// that describes the probability of `k` successes (random draws for which the object /// drawn has a specified feature) in `n` draws, without replacement, from a finite diff --git a/src/FSharp.Stats/Distributions/Discrete/NegativeBinomial.fs b/src/FSharp.Stats/Distributions/Discrete/NegativeBinomial.fs index 9b69b855..09986b46 100644 --- a/src/FSharp.Stats/Distributions/Discrete/NegativeBinomial.fs +++ b/src/FSharp.Stats/Distributions/Discrete/NegativeBinomial.fs @@ -374,11 +374,10 @@ type NegativeBinomial_failures = static member ToString r p = sprintf "NegativeBinomial_failures(r = %i, p = %f)" r p - /// Initializes a negative binomial distribution. - /// The negative binomial distribution is a discrete probability distribution
that models the number of failures needed k to get the rth success in repeated
independent Bernoulli trials with probability p.

The number of success states
The probability of each independent bernoulli trial
The number of failures before the rth success + /// Initializes a negative binomial distribution.
The negative binomial distribution is a discrete probability distribution
that models the number of failures needed k to get the rth success in repeated
independent Bernoulli trials with probability p.
 /// - /// - /// + /// The number of success states + /// The probability of each independent bernoulli trial /// /// /// diff --git a/src/FSharp.Stats/Fitting/CrossValidation.fs b/src/FSharp.Stats/Fitting/CrossValidation.fs index b7c6a59a..c321e432 100644 --- a/src/FSharp.Stats/Fitting/CrossValidation.fs +++ b/src/FSharp.Stats/Fitting/CrossValidation.fs @@ -153,6 +153,7 @@ module CrossValidation = /// Computes a repeated shuffel-and-split cross validation
p: percentage of training set size from original size,
iterations: number of random subset creation,
xData: rowwise x-coordinate matrix,
yData: yData vector
fit: x and y data lead to function that maps a xData row vector to a y-coordinate,
error: defines the error of the fitted y-coordinate and the actual y-coordinate,
getStDev: function that calculates the standard deviation from a seq<^T>. (Seq.stDev) /// /// + /// /// /// /// diff --git a/src/FSharp.Stats/Fitting/LinearRegression.fs b/src/FSharp.Stats/Fitting/LinearRegression.fs index 27be4c5c..66c50e4c 100644 --- a/src/FSharp.Stats/Fitting/LinearRegression.fs +++ b/src/FSharp.Stats/Fitting/LinearRegression.fs @@ -156,6 +156,7 @@ module LinearRegression = /// Returns the regression function of a line through the origin /// /// The functions slope + /// x value of which the corresponding y value should be predicted /// Function that takes a x value and returns the predicted y value /// /// diff --git a/src/FSharp.Stats/Fitting/LogisticRegression.fs b/src/FSharp.Stats/Fitting/LogisticRegression.fs index 95bf061f..296dabe8 100644 --- a/src/FSharp.Stats/Fitting/LogisticRegression.fs +++ b/src/FSharp.Stats/Fitting/LogisticRegression.fs @@ -190,6 +190,7 @@ module LogisticRegression = /// Returns the regression function /// /// + /// /// /// /// diff --git a/src/FSharp.Stats/ML/Imputation.fs b/src/FSharp.Stats/ML/Imputation.fs index 53ddae05..2ff7ea84 100644 --- a/src/FSharp.Stats/ML/Imputation.fs +++ b/src/FSharp.Stats/ML/Imputation.fs @@ -34,6 +34,8 @@ module Imputation = /// Imputation by random sampling from the input vector /// /// + /// + /// /// /// /// @@ -71,6 +73,9 @@ module Imputation = /// Imputation by k-nearest neighbour /// /// + /// + /// + /// /// /// /// diff --git a/src/FSharp.Stats/ML/Unsupervised/ClusterNumber.fs b/src/FSharp.Stats/ML/Unsupervised/ClusterNumber.fs index e9fe0579..e202526e 100644 --- a/src/FSharp.Stats/ML/Unsupervised/ClusterNumber.fs +++ b/src/FSharp.Stats/ML/Unsupervised/ClusterNumber.fs @@ -28,6 +28,8 @@ module ClusterNumber = /// Akaike Information Criterion (AIC) /// /// + /// + /// /// /// /// @@ -299,6 +301,7 @@ https://www.datanovia.com/en/lessons/determining-the-optimal-number-of-clusters- /// Generate uniform points within the range of `data`. /// /// + /// /// /// /// diff --git a/src/FSharp.Stats/ML/Unsupervised/KNN.fs b/src/FSharp.Stats/ML/Unsupervised/KNN.fs index 3711ba4d..83f39372 100644 --- a/src/FSharp.Stats/ML/Unsupervised/KNN.fs +++ b/src/FSharp.Stats/ML/Unsupervised/KNN.fs @@ -207,6 +207,7 @@ module KNN = /// Predict (or classify) the given point. /// /// the point to be classified. + /// /// /// /// // .. construct and fit the knnClassifier before .. @@ -221,6 +222,7 @@ module KNN = /// Predict (or classify) the given collection of points. /// /// the array of points to be classified. + /// /// /// /// // .. construct and fit the knnClassifier before .. diff --git a/src/FSharp.Stats/Testing/TestStatistics.fs b/src/FSharp.Stats/Testing/TestStatistics.fs index 2670a9fe..4c174f3d 100644 --- a/src/FSharp.Stats/Testing/TestStatistics.fs +++ b/src/FSharp.Stats/Testing/TestStatistics.fs @@ -10,15 +10,12 @@ module TestStatistics = open FSharp.Stats /// - /// Creates a new T-Test for a given statistic - /// with given degrees of freedom. + /// Represents the result of a T-Test with statistic, degrees of freedom, + /// and one-/two-tailed p-values. /// - /// - /// The test statistic. - /// The degrees of freedom for the numerator. - /// One Tailed/Sided. - /// One Tailed/Sided. - /// Two Tailed/Sided. + /// + /// Numeric type used for the test statistics. + /// type TTestStatistics<'T when 'T :> Numerics.INumber<'T> and Numerics.IFloatingPoint<'T> and Numerics.IExponentialFunctions<'T> @@ -27,10 +24,19 @@ module TestStatistics = and 'T: (static member op_Explicit: ^T -> float) and 'T : comparison > = { + /// The test statistic. Statistic : 'T + + /// The degrees of freedom for the numerator. DegreesOfFreedom : 'T + + /// One Tailed/Sided PValue. PValueLeft : 'T + + /// One Tailed/Sided PValue. PValueRight : 'T + + /// Two Tailed/Sided PValue. PValue : 'T } @@ -61,10 +67,10 @@ module TestStatistics = /// Creates a new F-Test for a given statistic /// with given degrees of freedom. /// - /// - /// The test statistic. - /// The degrees of freedom for the numerator. - /// The degrees of freedom for the denominator. + // /// + // /// The test statistic. + // /// The degrees of freedom for the numerator. + // /// The degrees of freedom for the denominator. type FTestStatistics<'T when 'T :> Numerics.INumber<'T> and Numerics.IFloatingPoint<'T> and Numerics.IExponentialFunctions<'T> @@ -73,9 +79,15 @@ module TestStatistics = and 'T: (static member op_Explicit: ^T -> float) and 'T : comparison > = { + /// The test statistic. Statistic : 'T + + /// The degrees of freedom for the numerator. DegreesOfFreedom1 : 'T + + /// The degrees of freedom for the denominator. DegreesOfFreedom2 : 'T + PValue : 'T PValueTwoTailed : 'T } @@ -107,12 +119,12 @@ module TestStatistics = /// Computes the Chi-Square test statistics for a given statistic /// with given degrees of freedom. /// - /// - /// The test statistic. - /// The degrees of freedom for the numerator. - /// One Tailed/Sided. - /// One Tailed/Sided. - /// Two Tailed/Sided. + // /// + // /// The test statistic. + // /// The degrees of freedom for the numerator. + // /// One Tailed/Sided. + // /// One Tailed/Sided. + // /// Two Tailed/Sided. type ChiSquareStatistics<'T when 'T :> Numerics.INumber<'T> and Numerics.IFloatingPoint<'T> and Numerics.IExponentialFunctions<'T> @@ -121,7 +133,9 @@ module TestStatistics = and 'T: (static member op_Explicit: ^T -> float) and 'T : comparison > = { + /// The test statistic. Statistic : 'T + /// The degrees of freedom for the numerator. DegreesOfFreedom : 'T /// one tailed/sided chiSquare pValue PValueLeft : 'T @@ -158,9 +172,9 @@ module TestStatistics = /// /// Computes the Wilcoxon test statistics for a given statistic. /// - /// The test statistic. - /// One Tailed/Sided. - /// Two Tailed/Sided. + // /// The test statistic. + // /// One Tailed/Sided. + // /// Two Tailed/Sided. type WilcoxonTestStatistics<'T when 'T :> Numerics.INumber<'T> and Numerics.IFloatingPoint<'T> and Numerics.IExponentialFunctions<'T> @@ -169,9 +183,13 @@ module TestStatistics = and 'T: (static member op_Explicit: ^T -> float) and 'T : comparison > = { + /// The test statistic. Statistic : 'T + /// One Tailed/Sided PValue. PValueLeft : 'T - PValueRight : 'T + /// One Tailed/Sided PValue. + PValueRight : 'T + /// Two Tailed/Sided PValue. PValueTwoTailed : 'T } let inline createWilcoxon<'T when 'T :> Numerics.INumber<'T>