Add Optimization

FSharp.Stats.sln

@@ -30,6 +30,7 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "content", "content", "{8E6D
 	ProjectSection(SolutionItems) = preProject
 		docs\content\BasicStats.fsx = docs\content\BasicStats.fsx
 		docs\content\Distributions.fsx = docs\content\Distributions.fsx
+		docs\content\Fitting.fsx = docs\content\Fitting.fsx
 		docs\content\index.fsx = docs\content\index.fsx
 		docs\content\Matrix_Vector.fsx = docs\content\Matrix_Vector.fsx
 		docs\content\Optimization.fsx = docs\content\Optimization.fsx

docs/content/BasicStats.fsx

@@ -14,6 +14,14 @@ open FSharp.Stats
 
 
 
+let v = 
+    vector [|2.0; 20.0; 1.|]
+
+
+Vector.interval v
+
+
+
 
 let inline stDevPopulationOfMean mean (items:seq<'T>) : 'U  =
     use e = items.GetEnumerator()
@@ -48,7 +56,7 @@ let inline stDevPopulation (items:seq<'T>) : 'U  =
 
 
 
-stDevPopulation [1.;2.]
+stDevPopulation [1.;2.;3.;4.;]
 
 let zero = LanguagePrimitives.GenericZero< float > 
 let one = LanguagePrimitives.GenericOne< float >
@@ -57,5 +65,3 @@ let one = LanguagePrimitives.GenericOne< float >
 
 LanguagePrimitives.GenericGreaterThan 2. 1.
 
-
-

docs/content/Fitting.fsx

@@ -0,0 +1,90 @@
+(*** hide ***)
+// This block of code is omitted in the generated HTML documentation. Use 
+// it to define helpers that you do not want to show in the documentation.
+#I "../../bin"
+#r "../../packages/build/FSharp.Plotly/lib/net40/Fsharp.Plotly.dll"
+open FSharp.Plotly
+(**
+Statistical testing
+===================
+
+FSharp.Stats a provides 
+
+
+
+*)
+#r "FSharp.Stats.dll"
+open FSharp.Stats
+open FSharp.Stats.Fitting
+
+(**
+Linear Regression
+-----------------
+*)
+
+
+
+
+// Test versus http://www.cyclismo.org/tutorial/R/linearLeastSquares.html
+let xVector = vector [2000.;   2001.;  2002.;  2003.;   2004.;]
+let yVector = vector [9.34;   8.50;  7.62;  6.93;  6.60;]
+
+let coeff   = Regression.Linear.coefficient xVector yVector
+let fit     = Regression.Linear.fit coeff
+let regLine = xVector |> Vector.map fit
+
+
+
+let summary = Regression.calulcateSumOfSquares fit xVector yVector
+
+let rsquared = Regression.calulcateDetermination summary
+
+let sigIntercept = Regression.ttestIntercept coeff.[0] summary
+let sigSlope     = Regression.ttestSlope coeff.[1] summary
+
+
+let anova = Regression.Linear.calculateANOVA coeff xVector yVector
+
+
+let aic = Regression.calcAIC 2. summary.Count summary.Error
+let bic = Regression.calcBIC 2. summary.Count summary.Error
+
+Regression.getResiduals fit xVector yVector
+Regression.calculateSSE fit xVector yVector
+
+(*** define-output:regression1 ***)
+[
+    Chart.Point(Seq.zip xVector yVector,Name="data points");
+    Chart.Line(Seq.zip xVector regLine,Name ="regression")
+]
+|> Chart.Combine
+(*** include-it:regression1 ***)
+
+
+
+
+
+
+
+let xVector' = vector [1290.;1350.;1470.;1600.;1710.;1840.;1980.;2230.;2400.;2930.;]
+let yVector' = vector [1182.;1172.;1264.;1493.;1571.;1711.;1804.;1840.;1956.;1954.;]
+
+
+let coeff'   = Regression.Polynomial.coefficient 2 xVector' yVector'
+
+let fit'     = Regression.Polynomial.fit 2 coeff'
+let regLine' = vector xVector' |> Vector.map fit'
+
+
+Regression.Polynomial.calculateANOVA 2 coeff' xVector' yVector'
+
+(*** define-output:polynomial1 ***)
+[
+    Chart.Point(Seq.zip xVector' yVector',Name="data points");
+    Chart.Spline(Seq.zip xVector' regLine',Name ="regression")
+]
+|> Chart.Combine
+(*** include-it:polynomial1 ***)
+
+
+

docs/content/Matrix_Vector.fsx

@@ -10,7 +10,7 @@ dfdf
 #r "FSharp.Stats.dll"
 //open FSharp
 open Microsoft.FSharp.Math
-
+open FSharp.Stats
 
 
 
@@ -37,6 +37,8 @@ let rv =
     rowvec [|2.0; 20.0; 1.|]
 
 
+
+
 //let f x = x**2.
 
 //"2.;3.;".Split(';') 

src/FSharp.Stats/Array.fs

@@ -0,0 +1,138 @@
+namespace FSharp.Stats
+
+
+/// Module to compute common statistical measure on array
+[<AutoOpen>]
+module Array =
+
+
+    let range (items:list<_>) =        
+        let rec loop l (minimum) (maximum) =
+            match l with
+            | h::t -> loop t (min h minimum) (max h maximum)
+            | [] -> Intervals.create minimum maximum          
+        //Init by fist value
+        match items with
+        | h::t  -> loop t h h
+        | [] -> Intervals.Interval.Empty
+
+
+    /// Arranges the items between the left and right border, that all items left of the pivot element are smaller and bigger on the right.
+    /// Function works in place and returns the index of the pivote element
+    let partionSortInPlace left right (items:array<'T>) =   
+
+        // Swaps items of left and right index
+        let swapInPlace left right (items:array<'T>) =
+            let tmp = items.[left]
+            items.[left]  <- items.[right]
+            items.[right] <- tmp
+
+        // Median of three optimization improves performance in general,
+        // and eliminates worst-case behavior for sorted or reverse-sorted data.    
+        let center = (right + left) / 2
+        if (items.[left] > items.[right]) then
+            swapInPlace left right items
+        if (items.[center] < items.[left]) then
+            swapInPlace center  left items
+        if (items.[center] > items.[right]) then
+            swapInPlace center right  items
+        // // pick the pivot point and save it
+        let pivot = items.[center]    
+
+        let rec moveRightIndex j =
+            if (pivot < items.[j]) then
+                moveRightIndex (j-1)
+            else
+                j
+
+        let rec moveLeftIndex i =
+            if (items.[i] < pivot) then
+                moveLeftIndex (i+1) 
+            else
+                i
+
+        let rec loop i j =
+            if (i < j) then
+                let j' = moveRightIndex j
+                let i' = moveLeftIndex i
+                if i' < j' then swapInPlace i' j' items
+
+                loop i' j'
+            else
+                i
+
+        let i = loop left right        
+        i
+
+
+    /// Finds the kth smallest element in an unordered array
+    let quickSelect k (items:array<'T>) =
+
+        let rec quickSelect' (items:array<'T>) left right k =
+
+            // get pivot position  
+            let pivot = partionSortInPlace left right items 
+
+            // if pivot is less than k, select from the right part  
+            if (pivot < k) then             
+                quickSelect' items (pivot + 1) right k
+            // if pivot is greater than k, select from the left side  
+            elif (pivot > k) then
+                quickSelect' items left (pivot - 1) k
+            // if equal, return the value
+            else 
+                items.[pivot]
+
+
+        let items' = Array.copy items
+        quickSelect' items' 0 (items'.Length - 1) (k - 1)
+
+
+    /// Computes the sample median
+    let inline median (items:array<'T>) =
+
+        // returns max element of array from index to right index
+        let rec max cMax index rigthIndex (input:array<'T>) =
+            if index <= rigthIndex then
+                let current = input.[index]
+                if cMax < current then
+                    max current (index+1) rigthIndex input
+                else
+                    max cMax (index+1) rigthIndex input
+            else
+                cMax
+
+
+        // Returns a tuple of two items whose mean is the median by quickSelect algorithm
+        let rec medianQuickSelect (items:array<'T>) left right k =
+
+            // get pivot position  
+            let pivotIndex = partionSortInPlace left right items 
+
+            // if pivot is less than k, select from the right part  
+            if (pivotIndex < k) then             
+                medianQuickSelect items (pivotIndex + 1) right k
+            // if pivot is greater than k, select from the left side  
+            elif (pivotIndex > k) then
+                medianQuickSelect items left (pivotIndex - 1) k
+            // if equal, return the value
+            else 
+                let n = items.Length
+                if n % 2 = 0 then
+                    (max (items.[pivotIndex-1]) 0 (pivotIndex-1) items,items.[pivotIndex])
+                else
+                    (items.[pivotIndex],items.[pivotIndex])
+
+
+        let zero = LanguagePrimitives.GenericZero< 'T > 
+        let one = LanguagePrimitives.GenericOne< 'T > 
+
+        if items.Length > 0 then
+            let items' = Array.copy items
+            let n = items'.Length
+            let mid = n / 2    
+            let m1,m2 = medianQuickSelect items' 0 (items'.Length - 1) (mid)
+            (m1 + m2) / (one + one)
+
+        else
+            zero / zero    

src/FSharp.Stats/Distributions/Continuous.fs

@@ -581,21 +581,21 @@ module Continuous =
 
 
     // F-distribution helper functions.
-    let fTCheckParam dof1 dof2 = if dof1 < 0.0 || dof2 < 0.0 then failwith "F-distribution should be parametrized by dof1 and dof2 > 0.0."
+    let fCheckParam dof1 dof2 = if dof1 < 0.0 || dof2 < 0.0 then failwith "F-distribution should be parametrized by dof1 and dof2 > 0.0."
 
     /// F-distribution
     type F =
         /// Computes the mean.
         static member Mean dof1 dof2 =
-            fTCheckParam dof1 dof2
+            fCheckParam dof1 dof2
             if dof2 <= 2. then
                 nan
             else
                 dof2 / (dof2 - 2.0)
 
         /// Computes the variance.
         static member Variance dof1 dof2 =
-            fTCheckParam dof1 dof2
+            fCheckParam dof1 dof2
             if dof2 <= 4. then
                 nan
             else
@@ -604,20 +604,20 @@ module Continuous =
 
         /// Computes the standard deviation.
         static member StandardDeviation dof1 dof2 =
-            fTCheckParam dof1 dof2
+            fCheckParam dof1 dof2
             sqrt (F.Variance dof1 dof2)
 
 
         /// Produces a random sample using the current random number generator (from GetSampleGenerator()).
         static member Sample dof1 dof2 =
-            fTCheckParam dof1 dof2
+            fCheckParam dof1 dof2
             let gamma1 = Gamma.Sample (dof1 / 2.0) 2.0
             let gamma2 = Gamma.Sample (dof2 / 2.0) 2.0
             gamma1 / gamma2
 
         /// Computes the probability density function.
         static member PDF dof1 dof2 x =
-            fTCheckParam dof1 dof2
+            fCheckParam dof1 dof2
             if (x <= 0.) then
                 0.
             else
@@ -627,7 +627,7 @@ module Continuous =
 
         /// Computes the cumulative distribution function.
         static member CDF dof1 dof2 x =
-            fTCheckParam dof1 dof2
+            fCheckParam dof1 dof2
             if (x <= 0.) then
                 1.
             else
@@ -645,7 +645,7 @@ module Continuous =
 
         /// Returns the support of the exponential distribution: (0., Positive Infinity).
         static member Support dof1 dof2 =
-            fTCheckParam dof1 dof2
+            fCheckParam dof1 dof2
             (0., System.Double.PositiveInfinity)
 
     /// Initializes a F-distribution         
@@ -661,6 +661,9 @@ module Continuous =
         }   
 
 
+
 // ######
 // ... distribution 
-// ######
+// ######
+
+