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..67be4665 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 @@ -154,7 +154,7 @@ let mlcm = [2; 0; 4] ] |> array2D - |> Matrix.Generic.ofArray2D + |> Matrix.ofArray2D ) ) (** 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..64911c18 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" @@ -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 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" diff --git a/src/FSharp.Stats/Distributions/Continuous/StudentT.fs b/src/FSharp.Stats/Distributions/Continuous/StudentT.fs index d3053739..dd5df8fb 100644 --- a/src/FSharp.Stats/Distributions/Continuous/StudentT.fs +++ b/src/FSharp.Stats/Distributions/Continuous/StudentT.fs @@ -18,9 +18,9 @@ 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 + 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. /// /// 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/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 @@ - + 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/SpecialFunctions/Beta.fs b/src/FSharp.Stats/SpecialFunctions/Beta.fs index ab83be60..4f88f0f2 100644 --- a/src/FSharp.Stats/SpecialFunctions/Beta.fs +++ b/src/FSharp.Stats/SpecialFunctions/Beta.fs @@ -2,14 +2,21 @@ open System open FSharp.Stats +open FsMath +open FsMath.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 = 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 + + 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) @@ -19,7 +26,19 @@ 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<'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 +48,19 @@ 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<'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 +71,19 @@ 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<'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 +94,19 @@ 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<'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 /// @@ -79,29 +134,29 @@ module Beta = let c = 1.0 let d = let tmp = 1.0 - (qab * x / qap) - if (abs tmp < FPMIN) then 1. / FPMIN else 1. / tmp + if (abs tmp < f_FPMIN) then 1. / f_FPMIN else 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 d' = let tmp = 1.0 + (aa*d) - if (abs tmp < FPMIN) then 1. / FPMIN else 1. / tmp + if (abs tmp < f_FPMIN) then 1. / f_FPMIN else 1. / tmp let c' = let tmp = 1.0 + (aa/c) - if (abs tmp < FPMIN) then FPMIN else tmp + 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 < FPMIN) then 1. / FPMIN else 1. / tmp + if (abs tmp < f_FPMIN) then 1. / f_FPMIN else 1. / tmp let c'' = let tmp = 1.0 + (aa'/c') - if (abs tmp < FPMIN) then FPMIN else tmp + if (abs tmp < f_FPMIN) then f_FPMIN else tmp let del = d''*c'' let h'' = h' * del - if abs (del - 1.0) <= EPS then + if abs (del - 1.0) <= f_EPS then if isSymmetryTransformation then 1.0 - (bt*h''/a) else bt*h''/a else if m < 140 then diff --git a/src/FSharp.Stats/SpecialFunctions/Gamma.fs b/src/FSharp.Stats/SpecialFunctions/Gamma.fs index 189af2aa..724c9738 100644 --- a/src/FSharp.Stats/SpecialFunctions/Gamma.fs +++ b/src/FSharp.Stats/SpecialFunctions/Gamma.fs @@ -1,6 +1,8 @@ namespace FSharp.Stats.SpecialFunctions open System +open FsMath +open FsMath.GenericMath /// Approximations for the gamma function and related functions. /// @@ -23,15 +25,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 +62,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 +107,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 +128,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 +325,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 +343,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 +470,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) diff --git a/src/FSharp.Stats/Testing/TestStatistics.fs b/src/FSharp.Stats/Testing/TestStatistics.fs index 955c125f..4c174f3d 100644 --- a/src/FSharp.Stats/Testing/TestStatistics.fs +++ b/src/FSharp.Stats/Testing/TestStatistics.fs @@ -1,92 +1,213 @@ namespace FSharp.Stats.Testing +open System +open FsMath +open FsMath.GenericMath +// TODO: Update specific distributions to support generic type 'T to avoid explicit float casting 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. /// - 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 - } - - let createTTest statistic dof = - let cdf = Distributions.Continuous.StudentT.CDF 0. 1. dof statistic - let pvalue = if statistic > 0. then 1. - cdf else cdf - {Statistic=statistic; DegreesOfFreedom=dof; PValueLeft=1. - pvalue; PValueRight=pvalue; PValue=pvalue*2.;} + /// + /// Numeric type used for the test statistics. + /// + type TTestStatistics<'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 > = + { + /// 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 + } + + 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= 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 - } - - let createFTest statistic dof1 dof2 = - let cdf = Distributions.Continuous.F.CDF dof1 dof2 statistic + // /// + // /// 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 > = + { + /// 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 + } + + 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 - } - - - 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.} + // /// + // /// 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 > = + { + /// The test statistic. + Statistic : 'T + /// The degrees of freedom for the numerator. + 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 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 > = + { + /// The test statistic. + Statistic : 'T + /// One Tailed/Sided PValue. + PValueLeft : 'T + /// One Tailed/Sided PValue. + PValueRight : 'T + /// Two Tailed/Sided PValue. + 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 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"