From f612db93a1302f3c66967e01d8c5a0c48c47bade Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Tue, 4 May 2010 05:52:21 +0800 Subject: [PATCH] Fixed a bug in the StudentT sampler. Reworked the unit tests for sampling methods to be more accurate. --- .../Distributions/Continuous/StudentT.cs | 4 +- src/Numerics/Version.cs | 6 +- .../CommonDistributionTests.cs | 83 +++++++++---------- .../Multivariate/NormalGammaTests.cs | 64 ++++++-------- 4 files changed, 67 insertions(+), 90 deletions(-) diff --git a/src/Numerics/Distributions/Continuous/StudentT.cs b/src/Numerics/Distributions/Continuous/StudentT.cs index 1a30c051..a2b72107 100644 --- a/src/Numerics/Distributions/Continuous/StudentT.cs +++ b/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; } } } diff --git a/src/Numerics/Version.cs b/src/Numerics/Version.cs index 2deb9659..f9721d52 100644 --- a/src/Numerics/Version.cs +++ b/src/Numerics/Version.cs @@ -22,9 +22,9 @@ // OTHER DEALINGS IN THE SOFTWARE. // /* 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: AssemblyFileVersion("2010.04.23.604")] \ No newline at end of file +[assembly: AssemblyVersion("2010.05.3.1269")] +[assembly: AssemblyFileVersion("2010.05.3.1269")] \ No newline at end of file diff --git a/src/UnitTests/DistributionTests/CommonDistributionTests.cs b/src/UnitTests/DistributionTests/CommonDistributionTests.cs index a267efd4..454f535b 100644 --- a/src/UnitTests/DistributionTests/CommonDistributionTests.cs +++ b/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; - // The accuracy of the histograms. - public static double SampleAccuracy = 0.01; - // The number of buckets to use to test against the cdf. + public static int NumberOfTestSamples = 3500000; + // The number of buckets in the histogram for the sampling function tests. 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 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(); } - - 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()); - } + VapnikChervonenkisTest(Error, ErrorProbability, samples, cd); } } @@ -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); - - 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()); - } + VapnikChervonenkisTest(Error, ErrorProbability, dd.Samples().Select(x => (double) x).Take(NumberOfTestSamples), dd); } 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); - 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()); - } + /// + /// 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". + /// 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 . + /// Note that for intervals on the real line the VC-dimension is 2. + /// + /// The error we are willing to tolerate. + /// The error probability we are willing to tolerate. + /// The samples to use for testing. + /// The distribution we are testing. + public static void VapnikChervonenkisTest(double epsilon, double delta, IEnumerable 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()); } } } diff --git a/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs b/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs index bbb6605e..5777e610 100644 --- a/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs +++ b/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); - for (int i = 0; i < CommonDistributionTests.NumberOfBuckets; i++) - { - var bucket = histogram[i]; - double empiricalProbability = bucket.Count / (double)CommonDistributionTests.NumberOfTestSamples; - double realProbability = precMarginal.CumulativeDistribution(bucket.UpperBound) - - precMarginal.CumulativeDistribution(bucket.LowerBound); - Assert.LessThan(Math.Abs(empiricalProbability - realProbability), CommonDistributionTests.SampleAccuracy, cd.ToString()); - } + CommonDistributionTests.VapnikChervonenkisTest( + CommonDistributionTests.Error, + CommonDistributionTests.ErrorProbability, + precs, + precMarginal); + + // Check the mean distribution. + CommonDistributionTests.VapnikChervonenkisTest( + CommonDistributionTests.Error, + CommonDistributionTests.ErrorProbability, + means, + meanMarginal); } /// @@ -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); - for (int i = 0; i < CommonDistributionTests.NumberOfBuckets; i++) - { - var bucket = histogram[i]; - double empiricalProbability = bucket.Count / (double)CommonDistributionTests.NumberOfTestSamples; - double realProbability = precMarginal.CumulativeDistribution(bucket.UpperBound) - - precMarginal.CumulativeDistribution(bucket.LowerBound); - Assert.LessThan(Math.Abs(empiricalProbability - realProbability), CommonDistributionTests.SampleAccuracy, cd.ToString()); - } + CommonDistributionTests.VapnikChervonenkisTest( + CommonDistributionTests.Error, + CommonDistributionTests.ErrorProbability, + precs, + precMarginal); + + // Check the mean distribution. + CommonDistributionTests.VapnikChervonenkisTest( + CommonDistributionTests.Error, + CommonDistributionTests.ErrorProbability, + means, + meanMarginal); } } } \ No newline at end of file