diff --git a/src/Numerics/Distributions/Discrete/Binomial.cs b/src/Numerics/Distributions/Discrete/Binomial.cs index 89d58544..720dfb35 100644 --- a/src/Numerics/Distributions/Discrete/Binomial.cs +++ b/src/Numerics/Distributions/Discrete/Binomial.cs @@ -255,8 +255,13 @@ namespace MathNet.Numerics.Distributions return 1.0; } - int k = (int) Math.Floor(x); - return (_n - k) * Combinatorics.Combinations(_n,k) * SpecialFunctions.BetaRegularized(_n - k, 1 + k, 1-_p); + double cdf = 0.0; + for (int i = 0; i <= (int) Math.Floor(x); i++) + { + cdf += Combinatorics.Combinations(_n, i) * Math.Pow(_p, i) * Math.Pow(1.0 - _p, _n - i); + } + + return cdf; } #endregion @@ -366,7 +371,7 @@ namespace MathNet.Numerics.Distributions /// /// Samples a Binomially distributed random variable. /// - /// The number of successful trials. + /// The number of successes in trials. public int Sample() { return DoSample(RandomSource, _p, _n); @@ -375,7 +380,7 @@ namespace MathNet.Numerics.Distributions /// /// Samples an array of Bernoulli distributed random variables. /// - /// a sequence of successful trial counts. + /// a sequence of successes in trials. public IEnumerable Samples() { while (true) @@ -392,7 +397,7 @@ namespace MathNet.Numerics.Distributions /// The random number generator to use. /// The success probability of a trial; must be in the interval [0.0, 1.0]. /// The number of trials; must be positive. - /// The number of successes in trials. + /// The number of successes in trials. public static int Sample(System.Random rnd, double p, int n) { if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n)) @@ -409,7 +414,7 @@ namespace MathNet.Numerics.Distributions /// The random number generator to use. /// The success probability of a trial; must be in the interval [0.0, 1.0]. /// The number of trials; must be positive. - /// a sequence of successful trial counts. + /// a sequence of successes in trials. public static IEnumerable Samples(System.Random rnd, double p, int n) { if (Control.CheckDistributionParameters && !IsValidParameterSet(p, n)) diff --git a/src/Numerics/Distributions/Discrete/Categorical.cs b/src/Numerics/Distributions/Discrete/Categorical.cs index d532e6cc..9cde7e6d 100644 --- a/src/Numerics/Distributions/Discrete/Categorical.cs +++ b/src/Numerics/Distributions/Discrete/Categorical.cs @@ -322,7 +322,7 @@ namespace MathNet.Numerics.Distributions /// The number of successful trials. public int Sample() { - return DoSample(RandomSource, _p); + return Sample(RandomSource, _p); } /// @@ -331,10 +331,7 @@ namespace MathNet.Numerics.Distributions /// a sequence of successful trial counts. public IEnumerable Samples() { - while (true) - { - yield return DoSample(RandomSource, _p); - } + return Samples(RandomSource, _p); } #endregion diff --git a/src/Numerics/Statistics/Histogram.cs b/src/Numerics/Statistics/Histogram.cs index 687d730c..c58e66c1 100644 --- a/src/Numerics/Statistics/Histogram.cs +++ b/src/Numerics/Statistics/Histogram.cs @@ -35,10 +35,10 @@ namespace MathNet.Numerics.Statistics /// /// A consists of a series of s, - /// each representing a region limited by a lower bound (inclusive) and an upper bound (exclusive). + /// each representing a region limited by a lower bound (exclusive) and an upper bound (inclusive). /// [Serializable] - public class Bucket : IComparable + public class Bucket : IComparable, ICloneable { /// /// This IComparer performs comparisons between a point and a bucket. @@ -55,11 +55,11 @@ namespace MathNet.Numerics.Statistics { if (bkt2.Width == 0.0) { - return -bkt1.Contains(bkt2.LowerBound); + return -bkt1.Contains(bkt2.UpperBound); } else { - return -bkt2.Contains(bkt1.LowerBound); + return -bkt2.Contains(bkt1.UpperBound); } } } @@ -108,6 +108,15 @@ namespace MathNet.Numerics.Statistics Count = count; } + /// + /// Creates a copy of the Bucket with the lowerbound, upperbound and counts exactly equal. + /// + /// A cloned Bucket object. + public Object Clone() + { + return new Bucket(LowerBound, UpperBound, Count); + } + /// /// Width of the Bucket. /// @@ -132,9 +141,9 @@ namespace MathNet.Numerics.Statistics /// smaller than the bucket, +1 if the point is larger than the bucket. public int Contains(double x) { - if (LowerBound <= x) + if (LowerBound < x) { - if (UpperBound > x) + if (UpperBound >= x) { return 0; } @@ -146,7 +155,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Comparison of two disjoint buckets. + /// Comparison of two disjoint buckets. The buckets cannot be overlapping. /// public int CompareTo(Bucket bucket) { @@ -200,7 +209,7 @@ namespace MathNet.Numerics.Statistics /// public override string ToString() { - return "[" + this.LowerBound + ";" + this.UpperBound + ")"; + return "(" + this.LowerBound + ";" + this.UpperBound + "] = " + this.Count; } } @@ -247,12 +256,13 @@ namespace MathNet.Numerics.Statistics double upper = data.Maximum(); double width = (upper - lower) / nbuckets; - // Add buckets for each bin; the biggest bucket must be slightly larger then the maximal element. - for (int n = 0; n < nbuckets-1; n++) + // Add buckets for each bin; the smallest bucket's lowerbound must be slightly smaller + // than the minimal element. + AddBucket(new Bucket(lower.Decrement(), lower + width)); + for (int n = 1; n < nbuckets; n++) { AddBucket(new Bucket(lower + n * width, lower + (n + 1) * width)); } - AddBucket(new Bucket(lower + (nbuckets - 1) * width, upper.Increment())); AddData(data); } @@ -294,20 +304,23 @@ namespace MathNet.Numerics.Statistics /// The datapoint which we want to add. public void AddData(double d) { + // Sort if needed. + LazySort(); + if (d < this.LowerBound) { - this[0].LowerBound = d; - this[0].Count++; + // Make the lower bound just slightly smaller than the datapoint so it is contained in this bucket. + buckets[0].LowerBound = d.Decrement(); + buckets[0].Count++; } else if (d > this.UpperBound) { - // 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++; + buckets[BucketCount - 1].UpperBound = d; + buckets[BucketCount - 1].Count++; } else { - GetBucketOf(d).Count++; + buckets[GetBucketIndexOf(d)].Count++; } } @@ -345,24 +358,14 @@ namespace MathNet.Numerics.Statistics } } - /// - /// Returns the Bucket with index . - /// - /// The index of the bucket to retrieve. - /// A reference to the bucket at index . - public Bucket GetBucket(int i) - { - return buckets[i]; - } - /// /// Returns the Bucket that contains the value v. /// /// The point to search the bucket for. - /// The bucket containing the point. + /// A copy of the bucket containing point . public Bucket GetBucketOf(double v) { - return buckets[GetBucketIndexOf(v)]; + return (Bucket) buckets[GetBucketIndexOf(v)].Clone(); } /// @@ -415,13 +418,13 @@ namespace MathNet.Numerics.Statistics /// Gets the n'th bucket. /// /// The index of the bucket to be returned. - /// The n'th bucket. + /// A copy of the n'th bucket. public Bucket this[int n] { get { LazySort(); - return buckets[n]; + return (Bucket) buckets[n].Clone(); } } diff --git a/src/UnitTests/DistributionTests/CommonDistributionTests.cs b/src/UnitTests/DistributionTests/CommonDistributionTests.cs index 2cc8e001..7b4abc92 100644 --- a/src/UnitTests/DistributionTests/CommonDistributionTests.cs +++ b/src/UnitTests/DistributionTests/CommonDistributionTests.cs @@ -32,8 +32,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests using System.Linq; using System.Collections.Generic; using MbUnit.Framework; - using MathNet.Numerics.Distributions; + using MathNet.Numerics.Random; using MathNet.Numerics.Statistics; + using MathNet.Numerics.Distributions; [TestFixture] public class CommonDistributionTests @@ -116,29 +117,34 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [MultipleAsserts] public void SampleFollowsCorrectDistribution() { + Random rnd = new MersenneTwister(); + // The test samples from the distributions, builds a histogram and checks // whether the histogram follows the CDF. foreach (var dd in discreteDistributions) { - int[] samples = new int[numberOfTestSamples]; + dd.RandomSource = rnd; + + double[] samples = new double[numberOfTestSamples]; for (int i = 0; i < numberOfTestSamples; i++) { - samples[i] = dd.Sample(); + samples[i] = (double) dd.Sample(); } - var histogram = new Histogram(samples.Select(x => (double)x), numberOfBuckets); + var histogram = new Histogram(samples, numberOfBuckets); for (int i = 0; i < numberOfBuckets; i++) { - var bucket = histogram.GetBucket(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); + Assert.LessThan(Math.Abs(empiricalProbability - realProbability), sampleAccuracy, dd.ToString()); } } foreach (var cd in continuousDistributions) { + cd.RandomSource = rnd; double[] samples = new double[numberOfTestSamples]; for (int i = 0; i < numberOfTestSamples; i++) { @@ -148,11 +154,11 @@ namespace MathNet.Numerics.UnitTests.DistributionTests var histogram = new Histogram(samples, numberOfBuckets); for (int i = 0; i < numberOfBuckets; i++) { - var bucket = histogram.GetBucket(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); + Assert.LessThan(Math.Abs(empiricalProbability - realProbability), sampleAccuracy, cd.ToString()); } } } @@ -161,35 +167,39 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [MultipleAsserts] public void SamplesFollowsCorrectDistribution() { + Random rnd = new MersenneTwister(); + // 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.GetBucket(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); + Assert.LessThan(Math.Abs(empiricalProbability - realProbability), sampleAccuracy, dd.ToString()); } } foreach (var cd in continuousDistributions) { + cd.RandomSource = rnd; var samples = cd.Samples().Take(numberOfTestSamples); var histogram = new Histogram(samples, numberOfBuckets); for (int i = 0; i < numberOfBuckets; i++) { - var bucket = histogram.GetBucket(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); + Assert.LessThan(Math.Abs(empiricalProbability - realProbability), sampleAccuracy, cd.ToString()); } } } diff --git a/src/UnitTests/StatisticsTests/HistogramTests.cs b/src/UnitTests/StatisticsTests/HistogramTests.cs index 5be89d2c..d06cdc22 100644 --- a/src/UnitTests/StatisticsTests/HistogramTests.cs +++ b/src/UnitTests/StatisticsTests/HistogramTests.cs @@ -95,7 +95,9 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests } [Test] + [Row(0.0, -1)] [Row(1.0, 0)] + [Row(1.5, 0)] [Row(2.0, 1)] [Row(-1.0, -1)] public void ValidateContains(double x, int r) @@ -106,9 +108,8 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests #endregion + #region Histogram Tests - - #region [Test] public void CanCreateEmptyHistogram() { @@ -123,9 +124,8 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests } [Test] - [Row(0.0, 0)] [Row(0.5, 0)] - [Row(1.0, 1)] + [Row(1.0, 0)] [Row(10.0, 3)] [Row(10000.0, 4)] public void CanGetBucketIndexOf(double x, double i) @@ -139,6 +139,21 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Assert.AreEqual(i, h.GetBucketIndexOf(x)); } + [Test] + [Row(0.0)] + [Row(-1.0)] + [ExpectedArgumentException] + public void CanGetBucketIndexOfFailsWhenBucketDoesntExist(double x) + { + var h = new Histogram(); + h.AddBucket(new Bucket(0.0, 1.0)); + h.AddBucket(new Bucket(1.0, 2.0)); + h.AddBucket(new Bucket(2.0, 3.0)); + h.AddBucket(new Bucket(3.0, 20.0)); + h.AddBucket(new Bucket(20.0, Double.PositiveInfinity)); + int i = h.GetBucketIndexOf(x); + } + [Test] public void CanGetBucketOf() { @@ -256,18 +271,19 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Console.WriteLine("{0}", hist); - for (int i = 0; i < 8; i++) + for (int i = 1; i < 9; 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(2.0, hist[0].Count); - Assert.AreEqual(0.5, hist.LowerBound); - Assert.AreEqual(9.5.Increment(), hist.UpperBound); + Assert.AreEqual(0.5.Decrement(), hist.LowerBound); + Assert.AreEqual(9.5, hist.UpperBound); } [Test] + [MultipleAsserts] public void SmallDatasetHistogramWithBounds() { Histogram hist = new Histogram(smallDataset, 10, 0.0, 10.0);