Add Generalised Linear Models to FSharpStats #334
Adress changes requested in #344
- Loading branch information
commit 8dcd8ab565b2a4ba53f0c6c718f62d7e4f5e1669
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -219,6 +219,96 @@ module LinearAlgebra = | |
| // else LinearAlgebraManaged.QR a | ||
| LinearAlgebraManaged.QR a | ||
|
|
||
| /// Performs QR decomposition using an alternative algorithm. | ||
| /// Returns the orthogonal matrix Q and the upper triangular matrix R. | ||
| 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 | ||
|
|
||
| /// Solves a linear system of equations using QR decomposition. | ||
| /// | ||
| /// Parameters: | ||
| /// - A: The coefficient matrix of the linear system. | ||
| /// - t: The target vector of the linear system. | ||
| /// | ||
| /// Returns: | ||
| /// - mX: The solution vector of the linear system. | ||
| /// - r: The upper triangular matrix obtained from QR decomposition. | ||
|
There was a problem hiding this comment. Thanks for the updates! FSharp.Stats/src/FSharp.Stats/Fitting/LinearRegression.fs Lines 128 to 145 in c5444cc The same for line 222 f. |
||
| 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), | ||
|
|
||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not sure if I get whats happening here:
iin line 262iwhich always should be 0, because it was initialized with 0. Just after the determination of qLength the value at this specific index is replaced in line 262.Obviously I miss a step, but is there a possibility to create the qLength by Vector.init (...) rather than keeping it mutable and accessing/mutating it at multiple positions? If not, I'm happy to merge it, but I got confused during my attempt to understand whats happening here..
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could you additionally set into -brackets a reference to this implementation and a short description what the benefit, drawback is when using this implementation over others.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I do not have a direct reference since I had a look at the mathematical explanation for QR via Gram Schmidt and tried implementing it myself. The only real upside using this method instead of the in FSharp.Stats established one is the difference in the output Dimensions in R. An actual reference would be the numpy implantation of the reduced qr https://numpy.org/doc/stable/reference/generated/numpy.linalg.qr.html .
I also think that this implemtation needs to be optimised in the future to perfect the GLM.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The mutation was done to simplify the function to match the mathematical explanation better. This is one of the things that could be updated when I get to it.