Add Generalised Linear Models to FSharpStats

docs/data/cheese.csv

@@ -0,0 +1,31 @@
+"","Taste","Acetic","H2S","Lactic"
+"1",12.3,94,23,0.86
+"2",20.9,174,155,1.53
+"3",39,214,230,1.57
+"4",47.9,317,1801,1.81
+"5",5.6,106,45,0.99
+"6",25.9,298,2000,1.09
+"7",37.3,362,6161,1.29
+"8",21.9,436,2881,1.78
+"9",18.1,134,47,1.29
+"10",21,189,65,1.58
+"11",34.9,311,465,1.68
+"12",57.2,630,2719,1.9
+"13",0.7,88,20,1.06
+"14",25.9,188,140,1.3
+"15",54.9,469,856,1.52
+"16",40.9,581,14589,1.74
+"17",15.9,120,50,1.16
+"18",6.4,224,110,1.49
+"19",18,190,480,1.63
+"20",38.9,230,8639,1.99
+"21",14,96,141,1.15
+"22",15.2,200,185,1.33
+"23",32,234,10322,1.44
+"24",56.7,349,26876,2.01
+"25",16.8,214,39,1.31
+"26",11.6,421,25,1.46
+"27",26.5,638,1056,1.72
+"28",0.7,206,50,1.25
+"29",13.4,331,800,1.08
+"30",5.5,481,120,1.25

src/FSharp.Stats/Algebra/LinearAlgebra.fs

@@ -219,6 +219,111 @@ module LinearAlgebra =
         //                    else LinearAlgebraManaged.QR a
         LinearAlgebraManaged.QR a
 
+    /// <summary>
+    /// Performs QR decomposition using an alternative algorithm.
+    /// QR decomposition is a method to decompose a matrix A into two components:
+    /// Q (an orthogonal matrix) and R (an upper triangular matrix),
+    /// such that A = Q * R. It is commonly used in solving linear systems,
+    /// least squares fitting, and eigenvalue problems.
+    /// </summary>
+    /// <returns>
+    /// A tuple containing:
+    /// <list type="bullet">
+    ///   <item><description>Q: The orthogonal matrix obtained from the decomposition.</description></item>
+    ///   <item><description>R: The upper triangular matrix obtained from the decomposition.</description></item>
+    /// </list>
+    /// </returns>
+    let qrAlternative (A: Matrix<float>) =
+        let m: int = A.NumRows
+        let n: int = A.NumCols
+
+        let q: Matrix<float> = Matrix.zero m n
+        let r: Matrix<float> = Matrix.zero n n
+        let qLengths: Vector<float> = Vector.zeroCreate n
+
+        let getVectorLength (v: Vector<float>) = Vector.fold (fun folder i -> folder+(i*i)) 0. v
+
+        let setqOfA (n: int) =
+            let aN: Vector<float> =  Matrix.getCol A n
+            let qN = 
+                if n = 0 then 
+                    aN 
+                else 
+                    Array.init (n) (fun i -> 
+                        let denominator = qLengths[i]
+                        let forNominator: Vector<float> = Matrix.getCol q i 
+                        let nominator: float = Vector.dot aN forNominator
+                        r.[i, n] <- nominator
+                        (nominator/denominator) * forNominator
+                    )
+                    |> Array.fold (fun folder  e -> folder-e ) aN
+            Matrix.setCol q n qN
+            qN  
+
+        for i=0 to n-1 do
+            let qN = setqOfA i 
+            let qLength = getVectorLength qN
+            let rValue = sqrt(qLength)
+            r[i,i] <- rValue
+            qLengths[i] <- qLength
+
+        for i=0 to n-1 do
+            let qN: Vector<float> = Matrix.getCol q i
+            let updateQ = (1./sqrt( qLengths[i]  )) * qN 
+            Matrix.setCol q i updateQ
+            for j=i+1 to n-1 do
+                let denominator = r[i, i]
+                let nominator = r[i, j]
+                r[i, j] <- (nominator/denominator)
+
+        q, r
+
+    /// <summary>
+    /// Solves a linear system of equations using QR decomposition.
+    /// </summary>
+    /// <param name="A">The coefficient matrix of the linear system.</param>
+    /// <param name="t">The target vector of the linear system.</param>
+    /// <returns>
+    /// A tuple containing:
+    /// <list type="bullet">
+    ///   <item><description>mX: The solution vector of the linear system.</description></item>
+    ///   <item><description>r: The upper triangular matrix obtained from QR decomposition.</description></item>
+    /// </list>
+    /// </returns>
+    let solveLinearQR (A: Matrix<float>) (t: Vector<float>) =
+        let m = A.NumRows
+        let n = A.NumCols
+
+        System.Diagnostics.Debug.Assert(m >= n) 
+
+        let q,r = qrAlternative A 
+
+        let QT = q.Transpose
+
+        let mX = Vector.zeroCreate n
+
+        let c: Vector<float> = QT * t
+
+        let rec build_mX_inner cross_prod i j =
+            if j=n then 
+                cross_prod
+            else
+                let newCrossprod = cross_prod + (r[i, j] * mX[j])
+                build_mX_inner newCrossprod i (j+1)
+
+        let rec build_mX_outer i =
+            if i<0 then 
+                ()
+            else
+                let crossProd = build_mX_inner 0. i (i+1)
+                mX[i] <- (c[i] - crossProd) / r[i, i]
+                build_mX_outer (i-1)
+
+        build_mX_outer (n-1)
+
+        mX,r
+
+
     ///Returns the full Singular Value Decomposition of the input MxN matrix 
     ///
     ///A : A = U * SIGMA * V**T in the tuple (S, U, V**T), 

src/FSharp.Stats/FSharp.Stats.fsproj

@@ -152,6 +152,7 @@
     <Compile Include="Fitting\LogisticRegression.fs" />
     <Compile Include="Fitting\QuantileNormalization.fs" />
     <Compile Include="Fitting\Spline.fs" />
+    <Compile Include="Fitting\GeneralisedLinearModel.fs" />
     <!-- ML -->
     <Compile Include="ML\SurprisalAnalysis.fs" />
     <Compile Include="ML\SimilarityMetrics.fs" />
@@ -170,7 +171,8 @@
     <None Include="Playground.fsx" />
   </ItemGroup>
   <ItemGroup>
-    <Content Include="../../lib/*.dll" PackagePath="netlib_LAPACK"></Content>
+    <Content Include="../../lib/*.dll" PackagePath="netlib_LAPACK">
+    </Content>
   </ItemGroup>
   <ItemGroup>
     <PackageReference Include="FSharpAux.Core" Version="2.0.0" />
@@ -182,4 +184,7 @@
   <ItemGroup>
     <None Include="../../README.md" Pack="true" PackagePath="\" />
   </ItemGroup>
-</Project>
+  <ItemGroup>
+    <InternalsVisibleTo Include="FSharp.Stats.Tests" />
+  </ItemGroup>
+</Project>