From 2d5561539b28a27d99703fd0f27c2834ab5e93f6 Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Sun, 23 Aug 2009 19:51:20 +0800 Subject: [PATCH] Added Multinomial distribution Added MersenneTwister Fixed namespace error in unit tests Signed-off-by: jvangael Signed-off-by: jvangael --- .../Distributions/Discrete/Bernoulli.cs | 2 +- .../Distributions/Multivariate/Dirichlet.cs | 6 +- .../Distributions/Multivariate/Multinomial.cs | 278 ++++++++++++++++ src/Numerics/Numerics.csproj | 2 + src/Numerics/Random/MersenneTwister.cs | 312 ++++++++++++++++++ .../Multivariate/MultinomialTests.cs | 117 +++++++ src/UnitTests/Random/MersenneTwisterTests.cs | 51 +++ src/UnitTests/Random/RandomTests.cs | 91 +++++ .../Random/SystemRandomExtensionTests.cs | 2 +- src/UnitTests/UnitTests.csproj | 3 + 10 files changed, 857 insertions(+), 7 deletions(-) create mode 100644 src/Numerics/Distributions/Multivariate/Multinomial.cs create mode 100644 src/Numerics/Random/MersenneTwister.cs create mode 100644 src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs create mode 100644 src/UnitTests/Random/MersenneTwisterTests.cs create mode 100644 src/UnitTests/Random/RandomTests.cs diff --git a/src/Numerics/Distributions/Discrete/Bernoulli.cs b/src/Numerics/Distributions/Discrete/Bernoulli.cs index 605c0294..43828368 100644 --- a/src/Numerics/Distributions/Discrete/Bernoulli.cs +++ b/src/Numerics/Distributions/Discrete/Bernoulli.cs @@ -54,7 +54,7 @@ namespace MathNet.Numerics.Distributions private Random _random; /// - /// Construct a new Bernoulli distribution. + /// Initializes a new instance of the Bernoulli class. /// /// The probability of generating one. /// If the Bernoulli parameter is not in the range [0,1]. diff --git a/src/Numerics/Distributions/Multivariate/Dirichlet.cs b/src/Numerics/Distributions/Multivariate/Dirichlet.cs index 8ecf7f88..df4f5af8 100644 --- a/src/Numerics/Distributions/Multivariate/Dirichlet.cs +++ b/src/Numerics/Distributions/Multivariate/Dirichlet.cs @@ -122,11 +122,7 @@ namespace MathNet.Numerics.Distributions throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } - _alpha = new double[alpha.Length]; - for (int i = 0; i < alpha.Length; i++) - { - _alpha[i] = alpha[i]; - } + _alpha = (double[]) alpha.Clone(); } /// diff --git a/src/Numerics/Distributions/Multivariate/Multinomial.cs b/src/Numerics/Distributions/Multivariate/Multinomial.cs new file mode 100644 index 00000000..0b3e00b3 --- /dev/null +++ b/src/Numerics/Distributions/Multivariate/Multinomial.cs @@ -0,0 +1,278 @@ +// +// 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. +// + +namespace MathNet.Numerics.Distributions +{ + using System; + using System.Collections.Generic; + using Properties; + + + /// + /// Implements the multinomial distribution. For details about this distribution, see + /// Wikipedia - Multinomial distribution. + /// + /// The distribution is parameterized by a vector of ratios: in other words, the parameter + /// does not have to be normalized and sum to 1. The reason is that some vectors can't be exactly normalized + /// to sum to 1 in floating point representation. + /// The distribution will use the by default. + /// Users can set the random number generator by using the property. + /// The statistics classes will check all the incoming parameters whether they are in the allowed + /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters + /// to false, all parameter checks can be turned off. + public class Multinomial + { + /// + /// Stores the normalized multinomial probabilities. + /// + private double[] _p; + + /// + /// The distribution's random number generator. + /// + private Random _random; + + /// + /// Initializes a new instance of the Multinomial class. + /// + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// If any of the probabilities are negative or do not sum to one. + public Multinomial(double[] p) + { + SetParameters(p); + RandomSource = new System.Random(); + } + + /* TODO + /// + /// Generate a multinomial distribution from histogram . The distribution will + /// not be automatically updated when the histogram changes. + /// + public Multinomial(Histogram h) + { + // The probability distribution vector. + _p = new double[h.BinCount]; + + // Fill in the distribution vector. + for (int i = 0; i < h.BinCount; i++) + { + _p[i] = h[i]; + } + + RandomNumberGenerator = new System.Random(); + }*/ + + /// + /// A string representation of the distribution. + /// + public override string ToString() + { + return "Multinomial(Dimension = " + _p.Length + ")"; + } + + /// + /// Checks whether the parameters of the distribution are valid. + /// + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// If any of the probabilities are negative returns false, or if the sum of parameters is 0.0; otherwise true + private static bool IsValidParameterSet(double[] p) + { + double sum = 0.0; + for (int i = 0; i < p.Length; i++) + { + if (p[i] < 0.0 || Double.IsNaN(p[i])) + { + return false; + } + else + { + sum += p[i]; + } + } + + if (sum == 0.0) + { + return false; + } + + return true; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// When the parameters don't pass the function. + private void SetParameters(double[] p) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _p = (double[])p.Clone(); + } + + /// + /// Gets or sets the probability of generating a one. + /// + public double[] P + { + get + { + return (double[]) _p.Clone(); + } + + set + { + SetParameters(value); + } + } + + /// + /// Gets or sets the random number generator which is used to draw random samples. + /// + public Random RandomSource + { + get + { + return _random; + } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + + _random = value; + } + } + + /// + /// Samples one multinomial distributed random variable; also known as the Discrete distribution. + /// + /// One random integer between 0 and the size of the multinomial (exclusive). + public int Sample() + { + return Sample(RandomSource, _p); + } + + /// + /// Samples a multinomially distributed random variable. + /// + /// The number of variables needed. + /// random integers between 0 and the size of the multinomial (exclusive). + public int[] Sample(int n) + { + return Sample(RandomSource, n, _p); + } + + /// + /// Samples one multinomial distributed random variable; also known as the Discrete distribution. + /// + /// The random number generator to use. + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// One random integer between 0 and the size of the multinomial (exclusive). + public static int Sample(System.Random rnd, double[] p) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + // The cumulative density of p. + double[] cp = UnnormalizedCDF(p); + + double u = rnd.NextDouble()*cp[cp.Length - 1]; + int idx = 0; + while (u > cp[idx]) + { + idx++; + } + return idx; + } + + /// + /// Samples a multinomially distributed random variable. + /// + /// The random number generator to use. + /// The number of variables needed. + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// random integers between 0 and the size of the multinomial (exclusive). + public static int[] Sample(System.Random rnd, int n, double[] p) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + // The cumulative density of p. + double[] cp = UnnormalizedCDF(p); + + int[] arr = new int[n]; + for (int i = 0; i < n; i++) + { + double u = rnd.NextDouble()*cp[cp.Length - 1]; + int idx = 0; + while (u > cp[idx]) + { + idx++; + } + arr[i] = idx; + } + + return arr; + } + + /// + /// Computes the unnormalized cumulative distribution function. This method performs no + /// parameter checking. + /// + /// An array of nonnegative ratios: this array does not need to be normalized + /// as this is often impossible using floating point arithmetic. + /// An array representing the unnormalized cumulative distribution function. + private static double[] UnnormalizedCDF(double[] p) + { + double[] cp = (double[]) p.Clone(); + + for (int i = 1; i < p.Length; i++) + { + cp[i] += cp[i - 1]; + } + + return cp; + } + } +} \ No newline at end of file diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index eb43cc7b..fee8bf91 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -61,6 +61,7 @@ + @@ -102,6 +103,7 @@ Resources.resx + diff --git a/src/Numerics/Random/MersenneTwister.cs b/src/Numerics/Random/MersenneTwister.cs new file mode 100644 index 00000000..c6a74cbc --- /dev/null +++ b/src/Numerics/Random/MersenneTwister.cs @@ -0,0 +1,312 @@ +// +// 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. +// + +/* + Original code's copyright and license: + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +namespace MathNet.Numerics.Random +{ + using System; + + /// + /// Random number generator using Mersenne Twister 19937 algorithm. + /// + public class MersenneTwister : AbstractRandomNumberGenerator, IDisposable + { + /// + /// Mersenne twister constant. + /// + private const uint _lower_mask = 0x7fffffff; + + /// + /// Mersenne twister constant. + /// + private const int _m = 397; + + /// + /// Mersenne twister constant. + /// + private const uint _matrix_a = 0x9908b0df; + + /// + /// Mersenne twister constant. + /// + private const int _n = 624; + + /// + /// Mersenne twister constant. + /// + private const double _reciprocal = 1.0/4294967295.0; + + /// + /// Mersenne twister constant. + /// + private const uint _upper_mask = 0x80000000; + + /// + /// Mersenne twister constant. + /// + private static readonly uint[] _mag01 = {0x0U, _matrix_a}; + + /// + /// Mersenne twister constant. + /// + private readonly uint[] _mt = new uint[624]; + + /// + /// Mersenne twister constant. + /// + private int mti = _n + 1; + + /// + /// Initializes a new instance of the class using + /// the current time as the seed. + /// + /// If the seed value is zero, it is set to one. Uses the + /// value of to + /// set whether the instance is thread safe. + public MersenneTwister() : this((int) DateTime.Now.Ticks) + { + } + + /// + /// Initializes a new instance of the class using + /// the current time as the seed. + /// + /// if set to true , the class is thread safe. + public MersenneTwister(bool threadSafe) : this((int) DateTime.Now.Ticks, threadSafe) + { + } + + /// + /// Initializes a new instance of the class. + /// + /// The seed value. + /// Uses the value of to + /// set whether the instance is thread safe. + public MersenneTwister(int seed) : this(seed, Control.ThreadSafeRandomNumberGenerators) + { + } + + /// + /// Initializes a new instance of the class. + /// + /// The seed value. + /// if set to true, the class is thread safe. + public MersenneTwister(int seed, bool threadSafe) : base(threadSafe) + { + init_genrand((uint)seed); + } + + /*/// + /// Initializes a new instance of the class. + /// + /// The initialization key. + public MersenneTwister(int[] init_key) + { + if (init_key == null) + { + throw new ArgumentNullException("init_key"); + } + uint[] array = new uint[init_key.Length]; + for (int i = 0; i < array.Length; i++) + { + array[i] = (uint) init_key[i]; + } + init_by_array(array); + } + */ + /* initializes _mt[_n] with a seed */ + + private void init_genrand(uint s) + { + _mt[0] = s & 0xffffffff; + for (mti = 1; mti < _n; mti++) + { + _mt[mti] = (1812433253*(_mt[mti - 1] ^ (_mt[mti - 1] >> 30)) + (uint) mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array _mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + _mt[mti] &= 0xffffffff; + /* for >32 bit machines */ + } + } + + + /* initialize by an array with array-length */ + /* init_key is the array for initializing keys */ + /* slight change for C++, 2004/2/26 */ + + /* private void init_by_array(uint[] init_key) + { + uint key_length = (uint) init_key.Length; + init_genrand(19650218); + uint i = 1; + uint j = 0; + uint k = (_n > key_length ? _n : key_length); + for (; k > 0; k--) + { + _mt[i] = (_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30))*1664525)) + init_key[j] + j; //non linear + _mt[i] &= 0xffffffff; // for WORDSIZE > 32 machines + i++; + j++; + if (i >= _n) + { + _mt[0] = _mt[_n - 1]; + i = 1; + } + if (j >= key_length) j = 0; + } + for (k = _n - 1; k > 0; k--) + { + _mt[i] = (_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30))*1566083941)) - i; // non linear + _mt[i] &= 0xffffffff; // for WORDSIZE > 32 machines + i++; + if (i >= _n) + { + _mt[0] = _mt[_n - 1]; + i = 1; + } + } + + _mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array + }*/ + + /* generates a random number on [0,0xffffffff]-interval */ + + private uint genrand_int32() + { + uint y; + + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + if (mti >= _n) + { + /* generate _n words at one time */ + int kk; + + if (mti == _n + 1) /* if init_genrand() has not been called, */ + init_genrand(5489); /* a default initial seed is used */ + + for (kk = 0; kk < _n - _m; kk++) + { + y = (_mt[kk] & _upper_mask) | (_mt[kk + 1] & _lower_mask); + _mt[kk] = _mt[kk + _m] ^ (y >> 1) ^ _mag01[y & 0x1]; + } + for (; kk < _n - 1; kk++) + { + y = (_mt[kk] & _upper_mask) | (_mt[kk + 1] & _lower_mask); + _mt[kk] = _mt[kk + (_m - _n)] ^ (y >> 1) ^ _mag01[y & 0x1]; + } + y = (_mt[_n - 1] & _upper_mask) | (_mt[0] & _lower_mask); + _mt[_n - 1] = _mt[_m - 1] ^ (y >> 1) ^ _mag01[y & 0x1]; + + mti = 0; + } + + y = _mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680; + y ^= (y << 15) & 0xefc60000; + y ^= (y >> 18); + + return y; + } + + /// + /// Returns a random number between 0.0 and 1.0. + /// + /// + /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. + /// + protected override double DoSample() + { + return genrand_int32() * _reciprocal; + } + + /* /// + /// Generates a random number on [0,1) with 53-bit resolution. + /// + /// A random number on [0,1) with 53-bit resolution. + public double NextDoubleResolution53() + { + ulong a = genrand_int32() >> 5, b = genrand_int32() >> 6; + return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0); + }*/ + + #region IDisposable Members + + /// + /// Performs application-defined tasks associated with freeing, releasing, or resetting unmanaged resources. + /// + public void Dispose() + { + //do nothing in the managed version. + } + + #endregion + } +} \ No newline at end of file diff --git a/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs b/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs new file mode 100644 index 00000000..7bfd0fdb --- /dev/null +++ b/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs @@ -0,0 +1,117 @@ +// +// 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. +// + +namespace MathNet.Numerics.UnitTests.DistributionTests +{ + using System; + using System.Linq; + using MbUnit.Framework; + using MathNet.Numerics.Distributions; + + [TestFixture] + public class MultinomialTests + { + double[] badP; + double[] badP2; + double[] smallP; + double[] largeP; + + [SetUp] + public void SetUp() + { + Control.CheckDistributionParameters = true; + badP = new double[] { -1.0, 1.0 }; + badP2 = new double[] { 0.0, 0.0 }; + smallP = new double[] {1.0, 1.0, 1.0}; + largeP = new double[] {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0}; + } + + [Test] + public void CanCreateMultinomial() + { + var m = new Multinomial(largeP); + AssertEx.AreEqual(largeP, m.P); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void MultinomialCreateFailsWithNegativeRatios() + { + var m = new Multinomial(badP); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void MultinomialCreateFailsWithAllZeroRatios() + { + var m = new Multinomial(badP2); + } + + [Test] + public void ValidateToString() + { + var b = new Multinomial(smallP); + AssertEx.AreEqual("Multinomial(Dimension = 3)", b.ToString()); + } + + [Test] + public void CanSetProbability() + { + var b = new Multinomial(largeP); + b.P = smallP; + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void SetProbabilityFails() + { + var b = new Multinomial(largeP); + b.P = badP; + } + + [Test] + public void CanSampleStatic() + { + var d = Multinomial.Sample(new Random(), largeP); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void FailSampleStatic() + { + var d = Multinomial.Sample(new Random(), badP); + } + + [Test] + public void CanSample() + { + var n = new Multinomial(largeP); + var d = n.Sample(); + } + } +} \ No newline at end of file diff --git a/src/UnitTests/Random/MersenneTwisterTests.cs b/src/UnitTests/Random/MersenneTwisterTests.cs new file mode 100644 index 00000000..6b5a4375 --- /dev/null +++ b/src/UnitTests/Random/MersenneTwisterTests.cs @@ -0,0 +1,51 @@ +// +// 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. +// + +namespace MathNet.Numerics.UnitTests.RandomTests +{ + using MbUnit.Framework; + using MathNet.Numerics.Random; + + [TestFixture] + public class MersenneTwisterTests : RandomTests + { + public MersenneTwisterTests() : base(typeof(MersenneTwister)) + { + } + + [Test, MultipleAsserts] + public void SampleKnownValues() + { + MersenneTwister mt = new MersenneTwister(0); + Assert.AreEqual(mt.NextDouble(), 0.5488135024320365); + Assert.AreEqual(mt.NextDouble(), 0.5928446165269344); + Assert.AreEqual(mt.NextDouble(), 0.7151893651381110); + Assert.AreEqual(mt.NextDouble(), 0.8442657442866512); + } + } +} \ No newline at end of file diff --git a/src/UnitTests/Random/RandomTests.cs b/src/UnitTests/Random/RandomTests.cs new file mode 100644 index 00000000..cde75317 --- /dev/null +++ b/src/UnitTests/Random/RandomTests.cs @@ -0,0 +1,91 @@ +// +// 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. +// + +namespace MathNet.Numerics.UnitTests.RandomTests +{ + using System; + using System.Threading; + using MbUnit.Framework; + using MathNet.Numerics.Random; + + public abstract class RandomTests + { + private const int _n = 10000; + private readonly Type _randomType; + + protected RandomTests(Type randomType) + { + _randomType = randomType; + } + + [Test] + public void Sample() + { + System.Random random = (System.Random) Activator.CreateInstance(_randomType, new object[] {false}); + double sum = 0; + for (int i = 0; i < _n; i++) + { + double next = random.NextDouble(); + sum += next; + Assert.IsTrue(next >= 0); + Assert.IsTrue(next <= 1); + } + //make sure are within 10% of the expected sum. + Assert.IsTrue(sum >= _n / 2.0 - .05 * _n); + Assert.IsTrue(sum <= _n / 2.0 + .05 * _n); + if (random is IDisposable) + { + ((IDisposable) random).Dispose(); + } + } + + [Test] + public void ThreadSafeSample() + { + System.Random random = (System.Random) Activator.CreateInstance(_randomType, new object[] {true}); + + Thread t1 = new Thread(runTest); + Thread t2 = new Thread(runTest); + t1.Start(random); + t2.Start(random); + t1.Join(); + t2.Join(); + } + + public void runTest(object random) + { + System.Random rng = (System.Random) random; + for (int i = 0; i < _n; i++) + { + double next = rng.NextDouble(); + Assert.IsTrue(next >= 0); + Assert.IsTrue(next <= 1); + } + } + } +} \ No newline at end of file diff --git a/src/UnitTests/Random/SystemRandomExtensionTests.cs b/src/UnitTests/Random/SystemRandomExtensionTests.cs index 20e066a3..8b76cefd 100644 --- a/src/UnitTests/Random/SystemRandomExtensionTests.cs +++ b/src/UnitTests/Random/SystemRandomExtensionTests.cs @@ -26,7 +26,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.UnitTests.DistributionTests +namespace MathNet.Numerics.UnitTests.RandomTests { using MbUnit.Framework; using MathNet.Numerics.Random; diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index f6cd4957..469019aa 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -72,6 +72,7 @@ + @@ -93,6 +94,8 @@ + +