Added Multinomial distribution Added MersenneTwister Fixed namespace error in unit tests

Signed-off-by: jvangael <jurgen.vangael@gmail.com> Signed-off-by: jvangael <jurgen.vangael@gmail.com>
17 years ago · 2d5561539b
10 changed files with 857 additions and 7 deletions
--- 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;

        /// <summary>
-        /// Construct a new Bernoulli distribution.
+        /// Initializes a new instance of the Bernoulli class.
        /// </summary>
        /// <param name="p">The probability of generating one.</param>
        /// <exception cref="ArgumentOutOfRangeException">If the Bernoulli parameter is not in the range [0,1].</exception>
--- 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();
        }

        /// <summary>
--- a/src/Numerics/Distributions/Multivariate/Multinomial.cs
+++ b/src/Numerics/Distributions/Multivariate/Multinomial.cs
@ -0,0 +1,278 @@
+// <copyright file="Bernoulli.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.Distributions
+{
+    using System;
+    using System.Collections.Generic;
+    using Properties;
+
+
+    /// <summary>
+    /// Implements the multinomial distribution. For details about this distribution, see 
+    /// <a href="http://en.wikipedia.org/wiki/Multinomial_distribution">Wikipedia - Multinomial distribution</a>.
+    /// </summary>
+    /// <remarks><para>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.</para>
+    /// <para>The distribution will use the <see cref="System.Random"/> by default. 
+    /// Users can set the random number generator by using the <see cref="RandomSource"/> property.</para>
+    /// <para>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.</para></remarks>
+    public class Multinomial
+    {
+        /// <summary>
+        /// Stores the normalized multinomial probabilities.
+        /// </summary>
+        private double[] _p;
+
+        /// <summary>
+        /// The distribution's random number generator.
+        /// </summary>
+        private Random _random;
+
+        /// <summary>
+        /// Initializes a new instance of the Multinomial class.
+        /// </summary>
+        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
+        /// as this is often impossible using floating point arithmetic.</param>
+        /// <exception cref="ArgumentException">If any of the probabilities are negative or do not sum to one.</exception>
+        public Multinomial(double[] p)
+        {
+            SetParameters(p);
+            RandomSource = new System.Random();
+        }
+
+        /* TODO
+        /// <summary>
+        /// Generate a multinomial distribution from histogram <paramref name="h"/>. The distribution will
+        /// not be automatically updated when the histogram changes.
+        /// </summary>
+        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();
+        }*/
+
+        /// <summary>
+        /// A string representation of the distribution.
+        /// </summary>
+        public override string ToString()
+        {
+            return "Multinomial(Dimension = " + _p.Length + ")";
+        }
+
+        /// <summary>
+        /// Checks whether the parameters of the distribution are valid. 
+        /// </summary>
+        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
+        /// as this is often impossible using floating point arithmetic.</param>
+        /// <returns>If any of the probabilities are negative returns false, or if the sum of parameters is 0.0; otherwise true</returns>
+        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;
+        }
+
+        /// <summary>
+        /// Sets the parameters of the distribution after checking their validity.
+        /// </summary>
+        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
+        /// as this is often impossible using floating point arithmetic.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the parameters don't pass the <see cref="IsValidParameterSet"/> function.</exception>
+        private void SetParameters(double[] p)
+        {
+            if (Control.CheckDistributionParameters && !IsValidParameterSet(p))
+            {
+                throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+            }
+
+            _p = (double[])p.Clone();
+        }
+
+        /// <summary>
+        /// Gets or sets the probability of generating a one.
+        /// </summary>
+        public double[] P
+        {
+            get
+            {
+                return (double[]) _p.Clone();
+            }
+
+            set
+            {
+                SetParameters(value);
+            }
+        }
+
+        /// <summary>
+        /// Gets or sets the random number generator which is used to draw random samples.
+        /// </summary>
+        public Random RandomSource
+        {
+            get
+            {
+                return _random;
+            }
+
+            set
+            {
+                if (value == null)
+                {
+                    throw new ArgumentNullException();
+                }
+
+                _random = value;
+            }
+        }
+
+        /// <summary>
+        /// Samples one multinomial distributed random variable; also known as the Discrete distribution.
+        /// </summary>
+        /// <returns>One random integer between 0 and the size of the multinomial (exclusive).</returns>
+        public int Sample()
+        {
+            return Sample(RandomSource, _p);
+        }
+
+        /// <summary>
+        /// Samples a multinomially distributed random variable.
+        /// </summary>
+        /// <param name="n">The number of variables needed.</param>
+        /// <returns><paramref name="n"/> random integers between 0 and the size of the multinomial (exclusive).</returns>
+        public int[] Sample(int n)
+        {
+            return Sample(RandomSource, n, _p);
+        }
+
+        /// <summary>
+        /// Samples one multinomial distributed random variable; also known as the Discrete distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
+        /// as this is often impossible using floating point arithmetic.</param>
+        /// <returns>One random integer between 0 and the size of the multinomial (exclusive).</returns>
+        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;
+        }
+
+        /// <summary>
+        /// Samples a multinomially distributed random variable.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="n">The number of variables needed.</param>
+        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
+        /// as this is often impossible using floating point arithmetic.</param>
+        /// <returns><paramref name="n"/> random integers between 0 and the size of the multinomial (exclusive).</returns>
+        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;
+        }
+
+        /// <summary>
+        /// Computes the unnormalized cumulative distribution function. This method performs no
+        /// parameter checking.
+        /// </summary>
+        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
+        /// as this is often impossible using floating point arithmetic.</param>
+        /// <returns>An array representing the unnormalized cumulative distribution function.</returns>
+        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;
+        }
+    }
+}
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@ -61,6 +61,7 @@
    <Compile Include="Distributions\IDiscreteDistribution.cs" />
    <Compile Include="Distributions\IDistribution.cs" />
    <Compile Include="Distributions\Multivariate\Dirichlet.cs" />
+    <Compile Include="Distributions\Multivariate\Multinomial.cs" />
    <Compile Include="IntegralTransforms\Algorithms\DiscreteHartleyTransform.Naive.cs" />
    <Compile Include="IntegralTransforms\Algorithms\DiscreteHartleyTransform.Options.cs" />
    <Compile Include="IntegralTransforms\HartleyOptions.cs" />
@ -102,6 +103,7 @@
      <DependentUpon>Resources.resx</DependentUpon>
    </Compile>
    <Compile Include="Random\AbstractRandomNumberGenerator.cs" />
+    <Compile Include="Random\MersenneTwister.cs" />
    <Compile Include="Random\SystemRandomExtensions.cs" />
    <Compile Include="Sampling\Sample.Random.cs" />
    <Compile Include="Sampling\Sample.Chebyshev.cs" />
--- a/src/Numerics/Random/MersenneTwister.cs
+++ b/src/Numerics/Random/MersenneTwister.cs
@ -0,0 +1,312 @@
+// <copyright file="MersenneTwister.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>
+
+/* 
+   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;
+
+    /// <summary>
+    /// Random number generator using Mersenne Twister 19937 algorithm.
+    /// </summary>
+    public class MersenneTwister : AbstractRandomNumberGenerator, IDisposable
+    {
+        /// <summary>
+        /// Mersenne twister constant.
+        /// </summary>
+        private const uint _lower_mask = 0x7fffffff;
+
+        /// <summary>
+        /// Mersenne twister constant.
+        /// </summary>
+        private const int _m = 397;
+
+        /// <summary>
+        /// Mersenne twister constant.
+        /// </summary>
+        private const uint _matrix_a = 0x9908b0df;
+
+        /// <summary>
+        /// Mersenne twister constant.
+        /// </summary>
+        private const int _n = 624;
+
+        /// <summary>
+        /// Mersenne twister constant.
+        /// </summary>
+        private const double _reciprocal = 1.0/4294967295.0;
+
+        /// <summary>
+        /// Mersenne twister constant.
+        /// </summary>
+        private const uint _upper_mask = 0x80000000;
+
+        /// <summary>
+        /// Mersenne twister constant.
+        /// </summary>
+        private static readonly uint[] _mag01 = {0x0U, _matrix_a};
+
+        /// <summary>
+        /// Mersenne twister constant.
+        /// </summary>
+        private readonly uint[] _mt = new uint[624];
+
+        /// <summary>
+        /// Mersenne twister constant.
+        /// </summary>
+        private int mti = _n + 1;
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="MersenneTwister"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        /// <remarks>If the seed value is zero, it is set to one. Uses the
+        /// value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
+        /// set whether the instance is thread safe.</remarks>
+        public MersenneTwister() : this((int) DateTime.Now.Ticks)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="MersenneTwister"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public MersenneTwister(bool threadSafe) : this((int) DateTime.Now.Ticks, threadSafe)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="MersenneTwister"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <remarks>Uses the value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
+        /// set whether the instance is thread safe.</remarks>        
+        public MersenneTwister(int seed) : this(seed, Control.ThreadSafeRandomNumberGenerators)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="MersenneTwister"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <param name="threadSafe">if set to <c>true</c>, the class is thread safe.</param>
+        public MersenneTwister(int seed, bool threadSafe) : base(threadSafe)
+        {
+            init_genrand((uint)seed);
+        }
+        
+        /*/// <summary>
+        /// Initializes a new instance of the <see cref="MersenneTwister"/> class.
+        /// </summary>
+        /// <param name="init_key">The initialization key.</param>
+        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;
+        }
+
+        /// <summary>
+        /// Returns a random number between 0.0 and 1.0.
+        /// </summary>
+        /// <returns>
+        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
+        /// </returns>
+        protected override double DoSample()
+        {
+            return genrand_int32() * _reciprocal;
+        }
+
+       /* /// <summary>
+        /// Generates a random number on [0,1) with 53-bit resolution.
+        /// </summary>
+        /// <returns>A random number on [0,1) with 53-bit resolution.</returns>
+        public double NextDoubleResolution53()
+        {
+            ulong a = genrand_int32() >> 5, b = genrand_int32() >> 6;
+            return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
+        }*/
+
+        #region IDisposable Members
+
+        /// <summary>
+        /// Performs application-defined tasks associated with freeing, releasing, or resetting unmanaged resources.
+        /// </summary>
+        public void Dispose()
+        {
+            //do nothing in the managed version.
+        }
+
+        #endregion
+    }
+}
--- a/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs
+++ b/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs
@ -0,0 +1,117 @@
+// <copyright file="MultinomialTests.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.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<double[]>(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<string>("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();
+        }
+    }
+}
--- a/src/UnitTests/Random/MersenneTwisterTests.cs
+++ b/src/UnitTests/Random/MersenneTwisterTests.cs
@ -0,0 +1,51 @@
+// <copyright file="SystemRandomExtensionTests.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.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);
+        }
+    }
+}
--- a/src/UnitTests/Random/RandomTests.cs
+++ b/src/UnitTests/Random/RandomTests.cs
@ -0,0 +1,91 @@
+// <copyright file="SystemRandomExtensionTests.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.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);
+            }
+        }
+    }
+}
--- a/src/UnitTests/Random/SystemRandomExtensionTests.cs
+++ b/src/UnitTests/Random/SystemRandomExtensionTests.cs
@ -26,7 +26,7 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>

-namespace MathNet.Numerics.UnitTests.DistributionTests
+namespace MathNet.Numerics.UnitTests.RandomTests
 {
    using MbUnit.Framework;
    using MathNet.Numerics.Random;
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@ -72,6 +72,7 @@
    <Compile Include="DistributionTests\Continuous\NormalTests.cs" />
    <Compile Include="DistributionTests\Discrete\BernoulliTests.cs" />
    <Compile Include="DistributionTests\Multivariate\DirichletTests.cs" />
+    <Compile Include="DistributionTests\Multivariate\MultinomialTests.cs" />
    <Compile Include="IntegralTransformsTests\HartleyTest.cs" />
    <Compile Include="IntegralTransformsTests\FourierTest.cs" />
    <Compile Include="IntegralTransformsTests\MatchingNaiveTransformTest.cs" />
@ -93,6 +94,8 @@
    <Compile Include="NumberTheoryTests\IntegerTheoryTest.cs" />
    <Compile Include="PrecisionTest.cs" />
    <Compile Include="Properties\AssemblyInfo.cs" />
+    <Compile Include="Random\MersenneTwisterTests.cs" />
+    <Compile Include="Random\RandomTests.cs" />
    <Compile Include="Random\SystemRandomExtensionTests.cs" />
    <Compile Include="SortingTests.cs" />
    <Compile Include="SpecialFunctionsTests\ErfTests.cs" />