Fixed a bug in the StudentT sampler.

Reworked the unit tests for sampling methods to be more accurate.

src/Numerics/Distributions/Continuous/StudentT.cs

@@ -473,9 +473,9 @@ namespace MathNet.Numerics.Distributions
         {
             double dummy = 0.0;
             var n = Normal.SampleBoxMuller(rnd, out dummy);
-            var g = Gamma.Sample(rnd, dof / 2.0, 1.0);
+            var g = Gamma.Sample(rnd, 0.5 * dof, 0.5);
 
-            return Math.Sqrt(2.0 * dof / g) * n;
+            return Math.Sqrt(dof / g) * n;
         }
     }
 }

src/Numerics/Version.cs

@@ -22,9 +22,9 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 /* This file is automatically generated - do not modify it. Change Version.tt instead.
-   Last generated on: 4/23/2010 10:04:41 AM
+   Last generated on: 03/05/2010 21:09:03
 */
 using System.Reflection;
 
-[assembly: AssemblyVersion("2010.04.23.604")]
+[assembly: AssemblyVersion("2010.05.3.1269")]
-[assembly: AssemblyFileVersion("2010.04.23.604")]
+[assembly: AssemblyFileVersion("2010.05.3.1269")]

src/UnitTests/DistributionTests/CommonDistributionTests.cs

@@ -44,11 +44,13 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
     public class CommonDistributionTests
     {
         // The number of samples we want.
-        public static int NumberOfTestSamples = 10000000;
+        public static int NumberOfTestSamples = 3500000;
-        // The accuracy of the histograms.
+        // The number of buckets in the histogram for the sampling function tests.
-        public static double SampleAccuracy = 0.01;
-        // The number of buckets to use to test against the cdf.
         public static int NumberOfBuckets = 100;
+        // The error we want to tolerate for sampling functions.
+        public static double Error = 0.01;
+        // The error probability we want to tolerate for sampling functions.
+        public static double ErrorProbability = 0.001;
         // The list of discrete distributions which we test.
         private List<IDiscreteDistribution> discreteDistributions;
         // The list of continuous distributions which we test.
@@ -75,7 +77,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
 
         [Test]
         [MultipleAsserts]
-        public void ValidateThatUnivariateDistributionsHaveRandomSource(int i)
+        public void ValidateThatUnivariateDistributionsHaveRandomSource()
         {
             foreach(var dd in discreteDistributions)
             {
@@ -90,7 +92,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
 
         [Test]
         [MultipleAsserts]
-        public void CanSetRandomSource(int i)
+        public void CanSetRandomSource()
         {
             foreach(var dd in discreteDistributions)
             {
@@ -105,7 +107,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
 
         [Test]
         [MultipleAsserts]
-        public void FailSetRandomSourceWithNullReference(int i)
+        public void FailSetRandomSourceWithNullReference()
         {
             foreach(var dd in discreteDistributions)
             {
@@ -127,18 +129,15 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         {
             Random rnd = new MersenneTwister(1);
 
-            // The test samples from the distributions, builds a histogram and checks
-            // whether the histogram follows the CDF.
             foreach (var dd in discreteDistributions)
             {
                 dd.RandomSource = rnd;
-
                 double[] samples = new double[NumberOfTestSamples];
                 for (int i = 0; i < NumberOfTestSamples; i++)
                 {
-                    samples[i] = (double) dd.Sample();
+                    samples[i] = (double)dd.Sample();
                 }
-
+                VapnikChervonenkisTest(Error, ErrorProbability, samples, dd);
             }
 
             foreach (var cd in continuousDistributions)
@@ -149,16 +148,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
                 {
                     samples[i] = cd.Sample();
                 }
-
+                VapnikChervonenkisTest(Error, ErrorProbability, samples, cd);
-                var histogram = new Histogram(samples, NumberOfBuckets);
-                for (int i = 0; i < NumberOfBuckets; i++)
-                {
-                    var bucket = histogram[i];
-                    double empiricalProbability = bucket.Count / (double)NumberOfTestSamples;
-                    double realProbability = cd.CumulativeDistribution(bucket.UpperBound)
-                        - cd.CumulativeDistribution(bucket.LowerBound);
-                    Assert.LessThan(Math.Abs(empiricalProbability - realProbability), SampleAccuracy, cd.ToString());
-                }
             }
         }
 
@@ -171,38 +161,41 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
         {
             Random rnd = new MersenneTwister(1);
 
-            // The test samples from the distributions, builds a histogram and checks
-            // whether the histogram follows the CDF.
             foreach (var dd in discreteDistributions)
             {
                 dd.RandomSource = rnd;
-                var samples = dd.Samples().Take(NumberOfTestSamples).Select(x => (double)x);
+                VapnikChervonenkisTest(Error, ErrorProbability, dd.Samples().Select(x => (double) x).Take(NumberOfTestSamples), dd);
-
-                var histogram = new Histogram(samples, NumberOfBuckets);
-                for (int i = 0; i < NumberOfBuckets; i++)
-                {
-                    var bucket = histogram[i];
-                    double empiricalProbability = bucket.Count / (double)NumberOfTestSamples;
-                    double realProbability = dd.CumulativeDistribution(bucket.UpperBound)
-                        - dd.CumulativeDistribution(bucket.LowerBound);
-                    Assert.LessThan(Math.Abs(empiricalProbability - realProbability), SampleAccuracy, dd.ToString());
-                }
             }
 
             foreach (var cd in continuousDistributions)
             {
                 cd.RandomSource = rnd;
-                var samples = cd.Samples().Take(NumberOfTestSamples);
+                VapnikChervonenkisTest(Error, ErrorProbability, cd.Samples().Take(NumberOfTestSamples), cd);
+            }
+        }
 
-                var histogram = new Histogram(samples, NumberOfBuckets);
+        /// <summary>
-                for (int i = 0; i < NumberOfBuckets; i++)
+        /// <para>Using VC-dimension, we can bound the probability of making an error when estimating empirical probability
-                {
+        /// distributions. We are using Theorem 2.41 in "All Of Nonparametric Statistics".
-                    var bucket = histogram[i];
+        /// http://books.google.com/books?id=MRFlzQfRg7UC&lpg=PP1&dq=all%20of%20nonparametric%20statistics&pg=PA22#v=onepage&q=%22shatter%20coe%EF%AC%83cients%20do%20not%22&f=false .</para>
-                    double empiricalProbability = bucket.Count / (double)NumberOfTestSamples;
+        /// <para>Note that for intervals on the real line the VC-dimension is 2.</para>
-                    double realProbability = cd.CumulativeDistribution(bucket.UpperBound)
+        /// </summary>
-                        - cd.CumulativeDistribution(bucket.LowerBound);
+        /// <param name="epsilon">The error we are willing to tolerate.</param>
-                    Assert.LessThan(Math.Abs(empiricalProbability - realProbability), SampleAccuracy, cd.ToString());
+        /// <param name="delta">The error probability we are willing to tolerate.</param>
-                }
+        /// <param name="s">The samples to use for testing.</param>
+        /// <param name="dist">The distribution we are testing.</param>
+        public static void VapnikChervonenkisTest(double epsilon, double delta, IEnumerable<double> s, IDistribution dist)
+        {
+            double N = (double) s.Count();
+            Assert.GreaterThan(N, Math.Ceiling(32.0 * Math.Log(16.0 / delta) / epsilon / epsilon));
+
+            var histogram = new Histogram(s, NumberOfBuckets);
+
+            for (int i = 0; i < NumberOfBuckets; i++)
+            {
+                double p = dist.CumulativeDistribution(histogram[i].UpperBound) - dist.CumulativeDistribution(histogram[i].LowerBound);
+                double pe = histogram[i].Count / N;
+                Assert.LessThan(Math.Abs(p - pe), epsilon, dist.ToString());
             }
         }
     }

src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs

@@ -224,27 +224,19 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
             var meanMarginal = cd.MeanMarginal();
             var precMarginal = cd.PrecisionMarginal();
 
-            // Check the mean distribution.
-            var histogram = new Histogram(means, CommonDistributionTests.NumberOfBuckets);
-            for (int i = 0; i < CommonDistributionTests.NumberOfBuckets; i++)
-            {
-                var bucket = histogram[i];
-                double empiricalProbability = bucket.Count / (double)CommonDistributionTests.NumberOfTestSamples;
-                double realProbability = meanMarginal.CumulativeDistribution(bucket.UpperBound)
-                    - meanMarginal.CumulativeDistribution(bucket.LowerBound);
-                Assert.LessThan(Math.Abs(empiricalProbability - realProbability), CommonDistributionTests.SampleAccuracy, cd.ToString());
-            }
-
             // Check the precision distribution.
-            histogram = new Histogram(precs, CommonDistributionTests.NumberOfBuckets);
+            CommonDistributionTests.VapnikChervonenkisTest(
-            for (int i = 0; i < CommonDistributionTests.NumberOfBuckets; i++)
+                CommonDistributionTests.Error,
-            {
+                CommonDistributionTests.ErrorProbability,
-                var bucket = histogram[i];
+                precs,
-                double empiricalProbability = bucket.Count / (double)CommonDistributionTests.NumberOfTestSamples;
+                precMarginal);
-                double realProbability = precMarginal.CumulativeDistribution(bucket.UpperBound)
+
-                    - precMarginal.CumulativeDistribution(bucket.LowerBound);
+            // Check the mean distribution.
-                Assert.LessThan(Math.Abs(empiricalProbability - realProbability), CommonDistributionTests.SampleAccuracy, cd.ToString());
+            CommonDistributionTests.VapnikChervonenkisTest(
-            }
+                CommonDistributionTests.Error,
+                CommonDistributionTests.ErrorProbability, 
+                means,
+                meanMarginal);
         }
 
         /// <summary>
@@ -265,27 +257,19 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
             var meanMarginal = cd.MeanMarginal();
             var precMarginal = cd.PrecisionMarginal();
 
-            // Check the mean distribution.
-            var histogram = new Histogram(means, CommonDistributionTests.NumberOfBuckets);
-            for (int i = 0; i < CommonDistributionTests.NumberOfBuckets; i++)
-            {
-                var bucket = histogram[i];
-                double empiricalProbability = bucket.Count / (double)CommonDistributionTests.NumberOfTestSamples;
-                double realProbability = meanMarginal.CumulativeDistribution(bucket.UpperBound)
-                    - meanMarginal.CumulativeDistribution(bucket.LowerBound);
-                Assert.LessThan(Math.Abs(empiricalProbability - realProbability), CommonDistributionTests.SampleAccuracy, cd.ToString());
-            }
-
             // Check the precision distribution.
-            histogram = new Histogram(precs, CommonDistributionTests.NumberOfBuckets);
+            CommonDistributionTests.VapnikChervonenkisTest(
-            for (int i = 0; i < CommonDistributionTests.NumberOfBuckets; i++)
+                CommonDistributionTests.Error,
-            {
+                CommonDistributionTests.ErrorProbability,
-                var bucket = histogram[i];
+                precs,
-                double empiricalProbability = bucket.Count / (double)CommonDistributionTests.NumberOfTestSamples;
+                precMarginal);
-                double realProbability = precMarginal.CumulativeDistribution(bucket.UpperBound)
+
-                    - precMarginal.CumulativeDistribution(bucket.LowerBound);
+            // Check the mean distribution.
-                Assert.LessThan(Math.Abs(empiricalProbability - realProbability), CommonDistributionTests.SampleAccuracy, cd.ToString());
+            CommonDistributionTests.VapnikChervonenkisTest(
-            }
+                CommonDistributionTests.Error,
+                CommonDistributionTests.ErrorProbability, 
+                means,
+                meanMarginal);
         }
     }
 }