From 8b73646e7823523b2ad04a263795de45e6ba5804 Mon Sep 17 00:00:00 2001
From: Jurgen Van Gael <jurgen.vangael@gmail.com>
Date: Sun, 20 Sep 2009 00:04:58 +0800
Subject: [PATCH] Added Histogram example in F#. Removed AlmostZero from
 incomplete Gamma function. Ported Pearson correlation computation. Added more
 histogram unit tests. Fixed a bug in the histogram code.

Signed-off-by: jvangael <jurgen.vangael@gmail.com>

Signed-off-by: jvangael <jurgen.vangael@gmail.com>
---
 src/FSharpExamples/FSharpExamples.fsproj      |  1 +
 src/FSharpExamples/Histogram.fs               | 43 +++++++++++
 src/Numerics/Numerics.csproj                  |  1 +
 src/Numerics/SpecialFunctions/Gamma.cs        |  6 +-
 src/Numerics/Statistics/Correlation.cs        | 76 +++++++++++++++++++
 src/Numerics/Statistics/Histogram.cs          | 12 +--
 .../StatisticsTests/CorrelationTests.cs       | 70 +++++++++++++++++
 .../StatisticsTests/HistogramTests.cs         | 38 ++++++++++
 src/UnitTests/UnitTests.csproj                |  1 +
 9 files changed, 240 insertions(+), 8 deletions(-)
 create mode 100644 src/FSharpExamples/Histogram.fs
 create mode 100644 src/Numerics/Statistics/Correlation.cs
 create mode 100644 src/UnitTests/StatisticsTests/CorrelationTests.cs

diff --git a/src/FSharpExamples/FSharpExamples.fsproj b/src/FSharpExamples/FSharpExamples.fsproj
index 19695e9d..29aae21a 100644
--- a/src/FSharpExamples/FSharpExamples.fsproj
+++ b/src/FSharpExamples/FSharpExamples.fsproj
@@ -52,6 +52,7 @@
   <ItemGroup>
     <Compile Include="DenseVector.fs" />
     <Compile Include="Apply.fs" />
+    <Compile Include="Histogram.fs" />
   </ItemGroup>
   <Import Project="$(MSBuildExtensionsPath)\FSharp\1.0\Microsoft.FSharp.Targets" />
   <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
diff --git a/src/FSharpExamples/Histogram.fs b/src/FSharpExamples/Histogram.fs
new file mode 100644
index 00000000..65c58083
--- /dev/null
+++ b/src/FSharpExamples/Histogram.fs
@@ -0,0 +1,43 @@
+// <copyright file="Histogram.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+open MathNet.Numerics.Statistics
+
+/// The number of buckets to use in our histogram.
+let B = 9
+
+/// Create a small dataset.
+let data = [| 0.5; 1.5; 2.5; 3.5; 4.5; 5.5; 6.5; 7.5; 8.5; 9.5 |]
+
+/// A histogram with 9 buckets for this dataset.
+let hist = new Histogram(data, B)
+
+// Print some histogram information.
+printfn "Histogram.ToString(): %O" hist
+for i in 0 .. B-1 do
+    printfn "Bucket %d contains %f datapoints." i hist.[i].Count
\ No newline at end of file
diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 45eefabf..9200e049 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -121,6 +121,7 @@
     <Compile Include="SpecialFunctions\Factorial.cs" />
     <Compile Include="SpecialFunctions\Gamma.cs" />
     <Compile Include="SpecialFunctions\Stability.cs" />
+    <Compile Include="Statistics\Correlation.cs" />
     <Compile Include="Statistics\DescriptiveStatistics.cs" />
     <Compile Include="Statistics\Histogram.cs" />
     <Compile Include="Statistics\Statistics.cs" />
diff --git a/src/Numerics/SpecialFunctions/Gamma.cs b/src/Numerics/SpecialFunctions/Gamma.cs
index 6f19a5ad..7ae265ac 100644
--- a/src/Numerics/SpecialFunctions/Gamma.cs
+++ b/src/Numerics/SpecialFunctions/Gamma.cs
@@ -290,9 +290,9 @@ namespace MathNet.Numerics
                 throw new ArgumentOutOfRangeException("a,x", Properties.Resources.ArgumentNotNegative);
             }
 
-            if (Precision.AlmostZero(a))
+            if (a.AlmostEqual(0.0))
             {
-                if (Precision.AlmostZero(x))
+                if (x.AlmostEqual(0.0))
                 {
                     // either 0 or 1, depending on the limit direction
                     return Double.NaN;
@@ -301,7 +301,7 @@ namespace MathNet.Numerics
                 return 1d;
             }
 
-            if (Precision.AlmostZero(x))
+            if (x.AlmostEqual(0.0))
             {
                 return 0d;
             }
diff --git a/src/Numerics/Statistics/Correlation.cs b/src/Numerics/Statistics/Correlation.cs
new file mode 100644
index 00000000..aea15752
--- /dev/null
+++ b/src/Numerics/Statistics/Correlation.cs
@@ -0,0 +1,76 @@
+// <copyright file="Correlation.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Statistics
+{
+    using System;
+    using System.Collections.Generic;
+    using Properties;
+
+    /// <summary>
+    /// A class with correlation measures between two datasets.
+    /// </summary>
+    public static class Correlation
+    {
+        /// <summary>
+        /// Computes the Pearson product-moment correlation coefficient.
+        /// </summary>
+        /// <param name="dataA">Sample data A.</param>
+        /// <param name="dataB">Sample data B.</param>
+        /// <returns>The Pearson product-moment correlation coefficient.</returns>
+        public static double Pearson(IEnumerable<double> dataA, IEnumerable<double> dataB)
+        {
+            int n = 0;
+            double r = 0.0;
+            double meanA = dataA.Mean();
+            double meanB = dataB.Mean();
+            double sdevA = dataA.StandardDeviation();
+            double sdevB = dataB.StandardDeviation();
+
+            IEnumerator<double> ieA = dataA.GetEnumerator();
+            IEnumerator<double> ieB = dataB.GetEnumerator();
+
+            while (ieA.MoveNext())
+            {
+                if (ieB.MoveNext() == false)
+                {
+                    throw new ArgumentOutOfRangeException("Datasets dataA and dataB need to have the same length.");
+                }
+
+                n++;
+                r += (ieA.Current - meanA) * (ieB.Current - meanB) / (sdevA * sdevB);
+            }
+            if (ieB.MoveNext() == true)
+            {
+                throw new ArgumentOutOfRangeException("Datasets dataA and dataB need to have the same length.");
+            }
+
+            return r / (n - 1);
+        }
+    }
+}
diff --git a/src/Numerics/Statistics/Histogram.cs b/src/Numerics/Statistics/Histogram.cs
index 553156d0..2fe05db2 100644
--- a/src/Numerics/Statistics/Histogram.cs
+++ b/src/Numerics/Statistics/Histogram.cs
@@ -200,7 +200,7 @@ namespace MathNet.Numerics.Statistics
         /// <returns></returns>
         public override string ToString()
         {
-            return "[" + this.LowerBound + ";" + this.UpperBound + "]";
+            return "[" + this.LowerBound + ";" + this.UpperBound + ")";
         }
     }
     
@@ -244,14 +244,15 @@ namespace MathNet.Numerics.Statistics
             }
 
             double lower = data.Minimum();
-            double upper = data.Maximum().Increment();
+            double upper = data.Maximum();
             double width = (upper - lower) / nbuckets;
 
-            // Add buckets for each bin.
-            for (int n = 0; n < nbuckets; n++)
+            // Add buckets for each bin; the biggest bucket must be slightly larger then the maximal element.
+            for (int n = 0; n < nbuckets-1; n++)
             {
                 AddBucket(new Bucket(lower + n * width, lower + (n + 1) * width));
             }
+            AddBucket(new Bucket(lower + (nbuckets - 1) * width, upper.Increment()));
 
             AddData(data);
         }
@@ -300,7 +301,8 @@ namespace MathNet.Numerics.Statistics
             }
             else if (d > this.UpperBound)
             {
-                this[BucketCount - 1].UpperBound = d;
+                // Make the upper bound just slightly larger then the datapoint so it is contained in this bucket.
+                this[BucketCount - 1].UpperBound = d.Increment();
                 this[BucketCount - 1].Count++;
             }
             else
diff --git a/src/UnitTests/StatisticsTests/CorrelationTests.cs b/src/UnitTests/StatisticsTests/CorrelationTests.cs
new file mode 100644
index 00000000..3188bf0e
--- /dev/null
+++ b/src/UnitTests/StatisticsTests/CorrelationTests.cs
@@ -0,0 +1,70 @@
+// <copyright file="CorrelationTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.StatisticsTests
+{
+    using System;
+    using System.Linq;
+    using System.Collections.Generic;
+    using MbUnit.Framework;
+    using MathNet.Numerics.Statistics;
+    using MathNet.Numerics.LinearAlgebra.Double;
+
+    [TestFixture]
+    public class CorrelationTests
+    {
+        private readonly IDictionary<string, StatTestData> mData = new Dictionary<string, StatTestData>();
+
+        public CorrelationTests()
+        {
+            StatTestData lottery = new StatTestData("./data/NIST/Lottery.dat");
+            mData.Add("lottery", lottery);
+            StatTestData lew = new StatTestData("./data/NIST/Lew.dat");
+            mData.Add("lew", lew);
+        }
+
+        [Test]
+        public void PearsonCorrelationTest()
+        {
+            var dataA = mData["lottery"].Data.Take(200);
+            var dataB = mData["lew"].Data.Take(200);
+
+            double corr = Correlation.Pearson(dataA, dataB);
+            AssertHelpers.AlmostEqual(corr, -0.029470861580726, 13);
+        }
+
+        [Test]
+        public void PearsonCorrelationTest_Fail()
+        {
+            var dataA = mData["lottery"].Data;
+            var dataB = mData["lew"].Data;
+
+            Assert.Throws<ArgumentOutOfRangeException>(() => Correlation.Pearson(dataA, dataB));
+        }
+    }
+}
diff --git a/src/UnitTests/StatisticsTests/HistogramTests.cs b/src/UnitTests/StatisticsTests/HistogramTests.cs
index f22980f5..5be89d2c 100644
--- a/src/UnitTests/StatisticsTests/HistogramTests.cs
+++ b/src/UnitTests/StatisticsTests/HistogramTests.cs
@@ -37,6 +37,8 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests
     [TestFixture]
     public class HistogramTests
     {
+        private double[] smallDataset = { 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5 };
+
         #region BucketTests
         [Test]
         public void CanCreateEmptyBucket()
@@ -245,6 +247,42 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests
             Assert.AreEqual(3, h[0].Count);
         }
 
+        [Test]
+        public void SmallDatasetHistogramWithoutBounds()
+        {
+            Histogram hist = new Histogram(smallDataset, 9);
+
+            Assert.AreEqual(9, hist.BucketCount);
+
+            Console.WriteLine("{0}", hist);
+
+            for (int i = 0; i < 8; i++)
+            {
+                Console.WriteLine("{0} : {1}", i, hist[i].Count);
+                Assert.AreEqual(1.0, hist[i].Count);
+            }
+            Assert.AreEqual(2.0, hist[8].Count);
+
+            Assert.AreEqual(0.5, hist.LowerBound);
+            Assert.AreEqual(9.5.Increment(), hist.UpperBound);
+        }
+
+        [Test]
+        public void SmallDatasetHistogramWithBounds()
+        {
+            Histogram hist = new Histogram(smallDataset, 10, 0.0, 10.0);
+
+            Assert.AreEqual(10, hist.BucketCount);
+
+            for (int i = 0; i < 10; i++)
+            {
+                Assert.AreEqual(1.0, hist[i].Count);
+            }
+
+            Assert.AreEqual(0.0, hist.LowerBound);
+            Assert.AreEqual(10.0, hist.UpperBound);
+        }
+
         #endregion
     }
 }
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index 25bc2450..6a899b88 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -109,6 +109,7 @@
     <Compile Include="SpecialFunctionsTests\ErfTests.cs" />
     <Compile Include="SpecialFunctionsTests\FactorialTest.cs" />
     <Compile Include="SpecialFunctionsTests\SpecialFunctionsTests.cs" />
+    <Compile Include="StatisticsTests\CorrelationTests.cs" />
     <Compile Include="StatisticsTests\HistogramTests.cs" />
     <Compile Include="StatisticsTests\DescriptiveStatisticsTests.cs" />
     <Compile Include="StatisticsTests\StatisticsTests.cs" />