From aab6c0186bccc13d9aac855cac7ec7647c1c87ef Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Tue, 22 Sep 2009 15:22:28 +0800 Subject: [PATCH] Initial version of VectorNormal distribution. --- .../Multivariate/VectorNormal.cs | 283 ++++++++++++++++++ src/Numerics/Numerics.csproj | 1 + .../CommonDistributionTests.cs | 4 +- .../Multivariate/VectorNormalTests.cs | 213 +++++++++++++ src/UnitTests/UnitTests.csproj | 1 + 5 files changed, 501 insertions(+), 1 deletion(-) create mode 100644 src/Numerics/Distributions/Multivariate/VectorNormal.cs create mode 100644 src/UnitTests/DistributionTests/Multivariate/VectorNormalTests.cs diff --git a/src/Numerics/Distributions/Multivariate/VectorNormal.cs b/src/Numerics/Distributions/Multivariate/VectorNormal.cs new file mode 100644 index 00000000..f9e44c91 --- /dev/null +++ b/src/Numerics/Distributions/Multivariate/VectorNormal.cs @@ -0,0 +1,283 @@ +// +// 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 Properties; + using MathNet.Numerics.LinearAlgebra.Double; + + /// + /// This class implements functionality for the multivariate normal distribution. This distribution + /// is parameterized by a mean vector and a covariance matrix. + /// + /// The distribution will use the by default. + /// Users can get/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 VectorNormal + { + + + + + + + // The Dirichlet distribution parameters. + private double[] _alpha; + + /// + /// The distribution's random number generator. + /// + private Random _random; + + /// + /// Initializes a new instance of the Dirichlet class. The distribution will + /// be initialized with the default random number generator. + /// + /// An array with the Dirichlet parameters. + public Dirichlet(double[] alpha) + { + SetParameters(alpha); + RandomSource = new Random(); + } + + /// + /// Constructs a new symmetric Dirichlet distribution. The distribution will + /// be initialized with the default random number generator. + /// + /// The value of each parameter of the Dirichlet distribution. + /// The dimension of the Dirichlet distribution. + public Dirichlet(double alpha, int k) + { + // Create a parameter structure. + double[] parm = new double[k]; + for (int i = 0; i < k; i++) + { + parm[i] = alpha; + } + + SetParameters(parm); + RandomSource = new Random(); + } + + /// + /// Checks whether the parameters of the distribution are valid: no parameter can be less than zero and + /// at least one parameter should be larger than zero. + /// + /// The parameters of the Dirichlet distribution. + /// True when the parameters are valid, false otherwise. + public static bool IsValidParameterSet(double[] alpha) + { + bool allzero = true; + + for (int i = 0; i < alpha.Length; i++) + { + if (alpha[i] < 0.0) + { + return false; + } + else if (alpha[i] > 0.0) + { + allzero = false; + } + } + + if (allzero) + { + return false; + } + + return true; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// The parameters of the Dirichlet distribution. + /// When the parameters don't pass the function. + private void SetParameters(double[] alpha) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(alpha)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _alpha = (double[]) alpha.Clone(); + } + + /// + /// A string representation of the distribution. + /// + public override string ToString() + { + return "Dirichlet(Dimension = " + this.Dimension + ")"; + } + + /// + /// Gets the dimension of the Dirichlet distribution. + /// + public int Dimension + { + get { return _alpha.Length; } + } + + /// + /// Gets or sets the parameters of the Dirichlet distribution. + /// + public double[] Alpha + { + get + { + return _alpha; + } + set + { + SetParameters(value); + } + } + + /// + /// The sum of the Dirichlet parameters. + /// + private double AlphaSum + { + get + { + double s = 0.0; + for (int i = 0; i < _alpha.Length; i++) + { + s += _alpha[i]; + } + return s; + } + } + + /// + /// Gets the mean of the Dirichlet distribution. + /// + public double[] Mean + { + get + { + double sum = AlphaSum; + double[] parm = new double[Dimension]; + for (int i = 0; i < Dimension; i++) + { + parm[i] = _alpha[i] / sum; + } + + return parm; + } + } + + /// + /// Gets the variance of the Dirichlet distribution. + /// + public double[] Variance + { + get + { + double s = this.AlphaSum; + double[] v = new double[_alpha.Length]; + for (int i = 0; i < _alpha.Length; i++) + { + v[i] = _alpha[i]*(s - _alpha[i])/(s*s*(s + 1.0)); + } + return v; + } + } + + /// + /// 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 a Dirichlet distributed random vector. + /// + public double[] Sample() + { + return Sample(RandomSource, _alpha); + } + + /// + /// Samples a Dirichlet distributed random vector. + /// + /// The random number generator to use. + /// The Dirichlet distribution parameter. + public static double[] Sample(System.Random rnd, double[] alpha) + { + if (Control.CheckDistributionParameters && ! IsValidParameterSet(alpha)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + int n = alpha.Length; + double[] gv = new double[n]; + double sum = 0.0; + for (int i = 0; i < n; i++) + { + if (alpha[i] == 0.0) + { + gv[i] = 0.0; + } + else + { + gv[i] = Gamma.Sample(rnd, alpha[i], 1.0); + } + + sum += gv[i]; + } + + for (int i = 0; i < n; i++) + { + gv[i] /= sum; + } + + return gv; + } + } +} \ No newline at end of file diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 52b1c4ba..7c093512 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -61,6 +61,7 @@ + diff --git a/src/UnitTests/DistributionTests/CommonDistributionTests.cs b/src/UnitTests/DistributionTests/CommonDistributionTests.cs index 3d5ab6ff..34cb6479 100644 --- a/src/UnitTests/DistributionTests/CommonDistributionTests.cs +++ b/src/UnitTests/DistributionTests/CommonDistributionTests.cs @@ -41,13 +41,15 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [SetUp] public void SetupDistributions() { - dists = new IDistribution[5]; + dists = new IDistribution[7]; dists[0] = new Beta(1.0, 1.0); dists[1] = new ContinuousUniform(0.0, 1.0); dists[2] = new Gamma(1.0, 1.0); dists[3] = new Normal(0.0, 1.0); dists[4] = new Bernoulli(0.6); + dists[5] = new Weibull(1.0, 1.0); + dists[6] = new DiscreteUniform(1, 10); } [Test] diff --git a/src/UnitTests/DistributionTests/Multivariate/VectorNormalTests.cs b/src/UnitTests/DistributionTests/Multivariate/VectorNormalTests.cs new file mode 100644 index 00000000..6371dcaa --- /dev/null +++ b/src/UnitTests/DistributionTests/Multivariate/VectorNormalTests.cs @@ -0,0 +1,213 @@ +// +// 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 VectorNormalTests + { + [SetUp] + public void SetUp() + { + Control.CheckDistributionParameters = true; + } + + //[Test] + //[ExpectedException(typeof(ArgumentOutOfRangeException))] + //public void NormalConstructorFail() + //{ + // Matrix cov = new DenseMatrix(new double[,] { { 1.0, 1.0 }, { -1.0, 2.0 } }); + // Vector mean = new DenseVector(new double[] { 5.0, 5.0 }); + + // // Build a new vector normal distribution. + // VectorNormal normal = new VectorNormal(mean, cov); + //} + + [Test] + public void StandardNormal() + { + VectorNormal normal = new VectorNormal(5); + + // Test the mean. + for (int i = 0; i < 5; i++) + { + Assert.AreEqual(0.0, normal.Mean[i]); + } + + // Test the covariance. + for (int i = 0; i < 5; i++) + { + for (int j = 0; j < 5; j++) + { + if (i == j) + { + Assert.AreEqual(1.0, normal.Covariance[i, j]); + } + else + { + Assert.AreEqual(0.0, normal.Covariance[i, j]); + } + } + } + + // Test the pdf. + Assert.AreEqual(0.010105326013812, normal.Density(new DenseVector(5, 0.0)), mAcceptableError); + Assert.AreEqual(8.294956719377678e-004, normal.Density(new DenseVector(5, 1.0)), mAcceptableError); + + // Test the mode. + for (int i = 0; i < 5; i++) + { + Assert.AreEqual(0.0, normal.Mode[i]); + } + + // Test the median. + for (int i = 0; i < 5; i++) + { + Assert.AreEqual(0.0, normal.Median[i]); + } + + // Test the entropy. + Assert.AreEqual(7.094692666023364, normal.Entropy, mAcceptableError); + } + + [Test] + public void NormalFromCovariance() + { + Matrix cov = new DenseMatrix(new double[,] { { 1.0, 0.9 }, { 0.9, 1.0 } }); + Vector mean = new DenseVector(new double[] { 5.0, 5.0 }); + + // Check that these are valid mean and covariances. + Assert.DoesNotThrow(() => VectorNormal.CheckParameters(mean, cov)); + + // Build a new vector normal distribution. + VectorNormal normal = new VectorNormal(mean, cov); + + // Test the mean. + Assert.AreEqual(5.0, normal.Mean[0]); + Assert.AreEqual(5.0, normal.Mean[1]); + + // Test the covariance. + Assert.AreEqual(1.0, normal.Covariance[0, 0]); + Assert.AreEqual(0.9, normal.Covariance[0, 1]); + Assert.AreEqual(0.9, normal.Covariance[1, 0]); + Assert.AreEqual(1.0, normal.Covariance[1, 1]); + + // Test the mode. + Assert.AreEqual(5.0, normal.Mode[0]); + Assert.AreEqual(5.0, normal.Mode[1]); + + // Test the median. + Assert.AreEqual(5.0, normal.Median[0]); + Assert.AreEqual(5.0, normal.Median[1]); + + // Test the entropy. + Assert.AreEqual(2.007511462998520, normal.Entropy, mAcceptableError); + + // Get the RNG. + System.Random rnd = normal.RandomNumberGenerator; + } + + + + + + + + + [Test] + public void HasRandomSource(int i) + { + VectorNormal d = new VectorNormal(0.3, 5); + Assert.IsNotNull(d.RandomSource); + } + + [Test] + public void CanSetRandomSource(int i) + { + VectorNormal d = new VectorNormal(0.3, 5); + d.RandomSource = new Random(); + } + + [Test] + [ExpectedException(typeof(ArgumentNullException))] + public void FailSetRandomSourceWithNullReference(int i) + { + VectorNormal d = new VectorNormal(0.3, 5); + d.RandomSource = null; + } + + [Test] + public void CanGetDimension() + { + VectorNormal d = new VectorNormal(0.3, 10); + Assert.AreEqual(10, d.Dimension); + } + + [Test] + public void ValidateMean() + { + VectorNormal d = new VectorNormal(0.3, 5); + + for (int i = 0; i < 5; i++) + { + AssertHelpers.AlmostEqual(0.3/1.5, d.Mean[i], 15); + } + } + + [Test] + public void ValidateVariance() + { + double[] alpha = new double[10]; + double sum = 0.0; + for (int i = 0; i < 10; i++) + { + alpha[i] = i; + sum += i; + } + + VectorNormal d = new VectorNormal(alpha); + + for (int i = 0; i < 10; i++) + { + AssertHelpers.AlmostEqual(i * (sum - i) / (sum * sum * (sum + 1.0)), d.Variance[i], 15); + } + } + + [Test] + public void CanSampleVectorNormal() + { + VectorNormal d = new VectorNormal(1.0, 5); + double[] s = d.Sample(); + } + } +} \ No newline at end of file diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 6a899b88..67a3b9cd 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -72,6 +72,7 @@ +