ADd ChiSquare contingency table

Author: muehlhaus

src/FSharp.Stats/AlgTypes/GenericMath.fs

@@ -61,7 +61,10 @@ module GenericMath =
     let inline floor<'T when 'T :> Numerics.IFloatingPoint<'T>> (x: 'T) : 'T =
         'T. Floor(x)
 
-
+    /// Generic floor function
+    let inline epsilon<'T when 'T :> Numerics.IFloatingPoint<'T>> () : 'T =
+        'T.CreateTruncating System.Double.Epsilon
+ 
     // let inline min x y = if x < y then x else y
     // let inline max x y = if x > y then x else y
 

src/FSharp.Stats/AlgTypes/SIMDUtils.fs

@@ -5,6 +5,9 @@ open System.Runtime.InteropServices
 open System.Runtime.CompilerServices
 
 
+
+
+
 /// <summary>
 /// A utility class for performing SIMD (Single Instruction, Multiple Data) operations on arrays.
 /// </summary>
@@ -155,6 +158,44 @@ type SIMDUtils() =
         results
 
 
+    static member inline mapFoldUnchecked<'T
+            when 'T :> Numerics.INumber<'T>
+            and 'T : (new: unit -> 'T)
+            and 'T : struct
+            and 'T :> ValueType>
+            (fMapSimd: Numerics.Vector<'T> -> Numerics.Vector<'T>,
+             fMap: 'T -> 'T,
+             fReduceSimd: Numerics.Vector<'T> -> Numerics.Vector<'T> -> Numerics.Vector<'T>,
+             fReduceScalar: 'T -> 'T -> 'T,
+             zero: 'T,
+             input: ReadOnlySpan<'T>)
+            : 'T =
+
+            let len = input.Length
+            let simdSize = Numerics.Vector<'T>.Count
+            let simdBlocks = len / simdSize
+            let tailStart = simdBlocks * simdSize
+
+            // SIMD map + accumulation
+            let mutable accVec = Numerics.Vector<'T>(zero)
+            let inputVec = MemoryMarshal.Cast<'T, Numerics.Vector<'T>>(input)
+
+            for i = 0 to simdBlocks - 1 do
+                let mapped = fMapSimd inputVec[i]
+                accVec <- fReduceSimd accVec mapped
+
+            // Reduce SIMD accumulator to scalar
+            let mutable acc = accVec.[0]
+            for i = 1 to simdSize - 1 do
+                acc <- fReduceScalar acc accVec.[i]
+
+            // Scalar tail
+            for i = tailStart to len - 1 do
+                acc <- fReduceScalar acc (fMap input[i])
+
+            acc
+
+
 
 type SIMDRangeUtils =
 

src/FSharp.Stats/AlgTypes/Vector.fs

@@ -211,15 +211,17 @@ type Vector =
             let slotCount = length / slotSize
             let ceiling = slotSize * slotCount
 
-            let mutable result = Numerics.Vector<'T>.Zero
             let v1Span = MemoryMarshal.Cast<'T, Numerics.Vector<'T>>(v1.AsSpan())
             let v2Span = MemoryMarshal.Cast<'T, Numerics.Vector<'T>>(v2.AsSpan())
-            let mutable scalarResult = LanguagePrimitives.GenericZero<'T>
+            let mutable result = Numerics.Vector<'T>.Zero
+            //let mutable scalarResult = LanguagePrimitives.GenericZero<'T>
 
             for i = 0 to slotCount - 1 do
                 result <- Numerics.Vector.Add(result, Numerics.Vector.Multiply(v1Span.[i], v2Span.[i])) 
-                scalarResult <- Numerics.Vector.Sum result
-        
+                //scalarResult <- Numerics.Vector.Sum result
+            
+            let mutable scalarResult = Numerics.Vector.Sum result
+
             for i = ceiling to v1.Length - 1 do
                 scalarResult <- scalarResult + (v1.[i] * v2.[i])
 

src/FSharp.Stats/Algebra/SVD.fs

@@ -1,6 +1,265 @@
 namespace FSharp.Stats.Algebra
 
+
+open System
 open FSharp.Stats
+open FSharp.Stats.Acceleration
+
+[<Struct>]
+type Householder<'T when 'T :> Numerics.INumber<'T>> =
+    {
+        V: Vector<'T>
+        Tau: 'T
+        Beta: 'T
+    }
+
+
+type Householder() =
+
+    static member inline create<'T when 'T :> Numerics.INumber<'T>
+                and 'T : (new: unit -> 'T)                              
+                and 'T : struct
+                and 'T : comparison
+                and 'T :> ValueType
+                and 'T :> Numerics.IRootFunctions<'T>>
+                (x: Vector<'T>) :  Householder<'T> =
+        
+        let xSpan = x.AsSpan()
+        let alpha = xSpan[0]
+        let tail = xSpan.Slice(1)
+
+        let zero = GenericMath.zero<'T>
+        let one = GenericMath.one<'T>
+
+        let sigma =
+            SIMDUtils.mapFoldUnchecked(
+                (fun v -> v * v),
+                (fun x -> x * x),
+                (+),
+                (+),
+                zero,
+                tail
+            )
+
+        if sigma.Equals(zero) then
+            let v = Vector.zeroCreate x.Length
+            v.[0] <- one
+            {
+                V = v
+                Tau = zero
+                Beta = alpha
+            }
+        else
+            let sum = alpha * alpha + sigma
+            let beta = GenericMath.sqrt sum
+
+            let v0 =
+                if alpha <= zero then alpha - beta
+                else -sigma / (alpha + beta)
+
+            let tau =
+                let v0Sq = v0 * v0
+                (v0Sq + v0Sq) / (sigma + v0Sq)
+
+            let v = Vector.divideScalar v0 x // SIMD-aware scalar division
+            v.[0] <- one
+
+            {
+                V = v
+                Tau = tau
+                Beta = beta
+            }
+
+    static member inline applyLeft<'T when 'T :> Numerics.INumber<'T>
+        and 'T : (new: unit -> 'T)
+        and 'T : struct
+        and 'T : equality
+        and 'T :> ValueType>
+        (h: Householder<'T>, A: Matrix<'T>, rowOffset: int) =
+
+        let v = h.V
+        let tau = h.Tau
+        let m = A.NumRows
+        let n = A.NumCols
+        let vLen = v.Length
+
+        for j = 0 to n - 1 do
+            let mutable dot = GenericMath.zero<'T>
+            for i = 0 to vLen - 1 do
+                dot <- dot + v.[i] * A.[rowOffset + i, j]
+            let scale = tau * dot
+            for i = 0 to vLen - 1 do
+                A.[rowOffset + i, j] <- A.[rowOffset + i, j] - scale * v.[i]
+
+
+    static member inline applyRight<'T when 'T :> Numerics.INumber<'T>
+        and 'T : (new: unit -> 'T)
+        and 'T : struct
+        and 'T : equality
+        and 'T :> ValueType>
+        (h: Householder<'T>, A: Matrix<'T>, colOffset: int) =
+
+        let v = h.V
+        let tau = h.Tau
+        let m = A.NumRows
+        let vLen = v.Length
+
+        for i = 0 to m - 1 do
+            let mutable dot = GenericMath.zero<'T>
+            for j = 0 to vLen - 1 do
+                dot <- dot + A.[i, colOffset + j] * v.[j]
+            let scale = tau * dot
+            for j = 0 to vLen - 1 do
+                A.[i, colOffset + j] <- A.[i, colOffset + j] - scale * v.[j]
+
+type Bidiagonalization() =
+
+    static member inline bidiagonalizeInPlace<'T
+        when 'T :> Numerics.INumber<'T>
+        and 'T : (new: unit -> 'T)
+        and 'T : struct
+        and 'T : comparison
+        and 'T :> ValueType
+        and 'T :> Numerics.IRootFunctions<'T>>
+        (A: Matrix<'T>) : unit =
+
+        let m = A.NumCols
+        let n = A.NumRows
+        let minMN = min m n
+
+        for k = 0 to minMN - 1 do
+            // --- LEFT REFLECTION: Column k (zero below diagonal) ---
+            let colLen = m - k
+            let colVector = Array.init colLen (fun i -> A.[k + i, k])
+            let hLeft = Householder.create colVector
+
+            Householder.applyLeft(hLeft, A, k)
+
+            // Overwrite A[k..,k] with Householder beta at top and zeros below
+            A.[k, k] <- hLeft.Beta
+            for i = k + 1 to m - 1 do
+                A.[i, k] <- GenericMath.zero<'T>
+
+            
+
+            // --- RIGHT REFLECTION: Row k (zero right of superdiagonal) ---
+            if k < n - 1 then
+                let rowLen = n - (k + 1)
+                let rowVector = Array.init rowLen (fun j -> A.[k, k + 1 + j])
+                let hRight = Householder.create rowVector
+
+                Householder.applyRight(hRight, A, k + 1)
+
+                // Overwrite A[k,k+1..] with Householder beta at front, zeros right
+                A.[k, k + 1] <- hRight.Beta
+                for j = k + 2 to n - 1 do
+                    A.[k, j] <- GenericMath.zero<'T>
+
+                
+[<Struct>]
+type Bidiagonal<'T when 'T :> Numerics.INumber<'T>> = {
+    D : Vector<'T>   // main diagonal
+    E : Vector<'T>   // superdiagonal (length n-1)
+}
+
+
+/// Givens Rotation (Generic) 
+module Givens =
+
+    let inline compute<'T
+        when 'T :> Numerics.INumber<'T>
+        and 'T : comparison
+        and 'T :> Numerics.IRootFunctions<'T>
+        and 'T :> Numerics.IFloatingPointIeee754<'T>>
+        (a: 'T) (b: 'T) : 'T * 'T =
+
+        if b = GenericMath.zero then GenericMath.one, GenericMath.zero
+        elif GenericMath.abs b > GenericMath.abs a then
+            let t = a / b
+            let s = GenericMath.one / GenericMath.sqrt(GenericMath.one + t * t)
+            s * t, s
+        else
+            let t = b / a
+            let c = GenericMath.one / GenericMath.sqrt(GenericMath.one + t * t)
+            c, c * t
+
+module GolubKahan =
+
+    let inline diagonalize<'T
+        when 'T :> Numerics.INumber<'T>
+        and 'T :> Numerics.IRootFunctions<'T>
+        and 'T :> Numerics.IFloatingPointIeee754<'T>
+        and 'T : comparison
+        and 'T : struct
+        and 'T : (new: unit -> 'T)
+        and 'T :> ValueType>
+        (b: Bidiagonal<'T>) : Vector<'T> =
+
+        let d = Array.copy b.D
+        let e = Array.copy b.E
+        let n = d.Length
+        let eps = GenericMath.epsilon()
+        let two = GenericMath.one + GenericMath.one
+
+        let mutable iter = 0
+        let maxIter = 1000
+        let mutable doneIterating = false
+
+       
+
+        while iter < maxIter && not doneIterating do
+            let mutable converged = true
+            
+            for i = 0 to n - 2 do
+                let tolerance = eps * (GenericMath.abs d.[i] + GenericMath.abs d.[i + 1])
+                if abs e.[i] > tolerance then
+                    converged <- false
+
+            if converged then
+                doneIterating <- true
+            else
+                // Wilkinson shift
+                let m = n - 1
+                let dm1 = d.[m - 1]
+                let dm  = d.[m]
+                let em1 = e.[m - 1]
+
+                let delta = (dm1 - dm) / two
+                let sign =
+                    if delta >= GenericMath.zero then GenericMath.one
+                    else -GenericMath.one
+
+                let denom = abs delta + sqrt (delta * delta + em1 * em1)
+                let mu = dm - sign * (em1 * em1) / denom
+
+                // Initial bulge
+                let mutable x = d.[0] * d.[0] - mu * mu
+                let mutable z = d.[0] * e.[0]
+
+                for k = 0 to n - 2 do
+                    let c, s = Givens.compute x z
+
+                    let dk  = d.[k]
+                    let ek  = e.[k]
+                    let dk1 = d.[k + 1]
+
+                    let tau1 = c * dk + s * ek
+                    let tau2 = -s * dk1
+
+                    d.[k]     <- c * tau1 + s * tau2
+                    e.[k]     <- c * ek - s * dk1
+                    d.[k + 1] <- s * tau1 - c * tau2
+
+                    if k < n - 2 then
+                        x <- e.[k]
+                        z <- -s * e.[k + 1]
+                        e.[k + 1] <- c * e.[k + 1]
+
+                iter <- iter + 1
+
+        d
+
+
 
 module SVD =
 

src/FSharp.Stats/FSharp.Stats.fsproj

@@ -56,6 +56,7 @@
     <Compile Include="Quantile.fs" />
     <Compile Include="Precision.fs" />
     <Compile Include="Geometry.fs" />
+    <Compile Include="Table.fs" />
     <!-- SpecialFunctions -->
     <Compile Include="SpecialFunctions\Gamma.fs" />
     <Compile Include="SpecialFunctions\Factorial.fs" />
@@ -123,6 +124,7 @@
     <!-- Testing -->
     <Compile Include="Testing\Tables.fs" />
     <Compile Include="Testing\TestStatistics.fs" />
+    <Compile Include="Testing\ContingencyTable.fs" />
     <Compile Include="Testing\ConfusionMatrix.fs" />
     <Compile Include="Testing\ComparisonMetrics.fs" />
     <Compile Include="Testing\Anova.fs" />