diff --git a/src/Numerics/Distributions/Discrete/Zipf.cs b/src/Numerics/Distributions/Discrete/Zipf.cs index 54232f75..34849ea7 100644 --- a/src/Numerics/Distributions/Discrete/Zipf.cs +++ b/src/Numerics/Distributions/Discrete/Zipf.cs @@ -238,14 +238,14 @@ namespace MathNet.Numerics.Distributions throw new NotSupportedException(); } - return (SpecialFunctions.GeneralHarmonic(_n, _s - 3) * Math.Pow(SpecialFunctions.GeneralHarmonic(_n, _s), 2) - SpecialFunctions.GeneralHarmonic(_n, _s - 1) * (3 * SpecialFunctions.GeneralHarmonic(_n, _s - 2) * SpecialFunctions.GeneralHarmonic(_n, _s) - Math.Pow(SpecialFunctions.GeneralHarmonic(_n, _s - 1), 2))) / Math.Pow(SpecialFunctions.GeneralHarmonic(_n, _s - 2) * SpecialFunctions.GeneralHarmonic(_n, _s) - Math.Pow(SpecialFunctions.GeneralHarmonic(_n, _s - 1), 2), 1.5); + return ((SpecialFunctions.GeneralHarmonic(_n, _s - 3) * Math.Pow(SpecialFunctions.GeneralHarmonic(_n, _s), 2)) - (SpecialFunctions.GeneralHarmonic(_n, _s - 1) * ((3 * SpecialFunctions.GeneralHarmonic(_n, _s - 2) * SpecialFunctions.GeneralHarmonic(_n, _s)) - Math.Pow(SpecialFunctions.GeneralHarmonic(_n, _s - 1), 2)))) / Math.Pow((SpecialFunctions.GeneralHarmonic(_n, _s - 2) * SpecialFunctions.GeneralHarmonic(_n, _s)) - Math.Pow(SpecialFunctions.GeneralHarmonic(_n, _s - 1), 2), 1.5); } } - /// + /// /// Computes the cumulative distribution function of the distribution. /// - /// The location at which to compute the cumulative density. + /// The integer location at which to compute the cumulative density. /// the cumulative density at . public double CumulativeDistribution(double x) { @@ -254,14 +254,7 @@ namespace MathNet.Numerics.Distributions return 0.0; } - var sum = 0.0; - var k = (int)Math.Ceiling(x) - 1; - for (var i = 0; i < k; i++) - { - sum += Math.Pow(i + 1, -_s); - } - - return sum / SpecialFunctions.GeneralHarmonic(_n, _s); + return SpecialFunctions.GeneralHarmonic((int)x, _s) / SpecialFunctions.GeneralHarmonic(_n, _s); } #endregion @@ -286,7 +279,7 @@ namespace MathNet.Numerics.Distributions { get { - throw new NotImplementedException(); + throw new NotSupportedException(); } } @@ -368,20 +361,25 @@ namespace MathNet.Numerics.Distributions /// a random number from the Zipf distribution. private static int DoSample(Random rnd, double s, int n) { - var r = rnd.NextDouble(); - var x = 1; + var r = 0.0; + while (r == 0.0) + { + r = rnd.NextDouble(); + } + var p = 1.0 / SpecialFunctions.GeneralHarmonic(n, s); - do + int i; + var sum = 0.0; + for (i = 1; i <= n; i++) { - r -= p; - if (r >= 0) + sum += p / Math.Pow(i, s); + if (sum >= r) { - x++; - p = p / Math.Pow((double)(x - 1) / x, s - 1); + break; } } - while (r >= 0); - return x; + + return i; } } } diff --git a/src/Numerics/Distributions/Multivariate/Dirichlet.cs b/src/Numerics/Distributions/Multivariate/Dirichlet.cs index 7b468e45..9e521685 100644 --- a/src/Numerics/Distributions/Multivariate/Dirichlet.cs +++ b/src/Numerics/Distributions/Multivariate/Dirichlet.cs @@ -215,6 +215,80 @@ namespace MathNet.Numerics.Distributions } } + /// + /// Gets the entropy of the distribution. + /// + public double Entropy + { + get + { + var num = _alpha.Sum(t => (t - 1) * SpecialFunctions.DiGamma(t)); + return SpecialFunctions.GammaLn(AlphaSum) + ((AlphaSum - Dimension) * SpecialFunctions.DiGamma(AlphaSum)) - num; + } + } + + /// + /// Computes the density of the distribution. + /// + /// The locations at which to compute the density. + /// the density at . + /// The Dirichlet distribution requires that the sum of the components of x equals 1. + /// You can also leave out the last component, and it will be computed from the others. + public double Density(double[] x) + { + return Math.Exp(DensityLn(x)); + } + + /// + /// Computes the log density of the distribution. + /// + /// The locations at which to compute the density. + /// the density at . + public double DensityLn(double[] x) + { + if (x == null) + { + throw new ArgumentNullException("x"); + } + + var flag = x.Length == (_alpha.Length - 1); + if ((x.Length != _alpha.Length) && !flag) + { + throw new ArgumentException("x"); + } + + var num = 0.0; + var num2 = 0.0; + for (var i = 0; i < x.Length; i++) + { + var d = x[i]; + if ((d <= 0.0) || (d >= 1.0)) + { + return 0.0; + } + + num += (_alpha[i] - 1.0) * Math.Log(d); + num2 += d; + } + + // Calculate x[Length - 1] element, if needed + if (flag) + { + if (num2 >= 1.0) + { + return 0.0; + } + + num += (_alpha[_alpha.Length - 1] - 1.0) * Math.Log(1.0 - num2); + } + else if (!num2.AlmostEqualInDecimalPlaces(1.0, 8)) + { + return 0.0; + } + + return -SpecialFunctions.GammaLn(AlphaSum) + num; + } + /// /// Gets or sets the random number generator which is used to draw random samples. /// diff --git a/src/Numerics/Distributions/Multivariate/InverseWishart.cs b/src/Numerics/Distributions/Multivariate/InverseWishart.cs index 2fa00128..6056adf3 100644 --- a/src/Numerics/Distributions/Multivariate/InverseWishart.cs +++ b/src/Numerics/Distributions/Multivariate/InverseWishart.cs @@ -200,6 +200,43 @@ namespace MathNet.Numerics.Distributions } } + /// + /// Gets the mode of the distribution. + /// + /// The mode of the distribution. + /// A. O'Hagan, and J. J. Forster (2004). Kendall's Advanced Theory of Statistics: Bayesian Inference. 2B (2 ed.). Arnold. ISBN 0-340-80752-0. + public Matrix Mode + { + get + { + return _s * (1.0 / (_nu + _s.RowCount + 1.0)); + } + } + + /// + /// Gets the variance of the distribution. + /// + /// The variance of the distribution. + /// Kanti V. Mardia, J. T. Kent and J. M. Bibby (1979). Multivariate Analysis. + public Matrix Variance + { + get + { + var res = _s.CreateMatrix(_s.RowCount, _s.ColumnCount); + for (var i = 0; i < res.RowCount; i++) + { + for (var j = 0; j < res.ColumnCount; j++) + { + var num1 = ((_nu - _s.RowCount + 1) * _s.At(i, j) * _s.At(i, j)) + ((_nu - _s.RowCount - 1) * _s.At(i, i) * _s.At(j, j)); + var num2 = (_nu - _s.RowCount) * (_nu - _s.RowCount - 1) * (_nu - _s.RowCount - 1) * (_nu - _s.RowCount - 3); + res.At(i, j, num1 / num2); + } + } + + return res; + } + } + /// /// Evaluates the probability density function for the inverse Wishart distribution. /// diff --git a/src/Numerics/Distributions/Multivariate/Multinomial.cs b/src/Numerics/Distributions/Multivariate/Multinomial.cs index d13acf0a..ba424745 100644 --- a/src/Numerics/Distributions/Multivariate/Multinomial.cs +++ b/src/Numerics/Distributions/Multivariate/Multinomial.cs @@ -28,6 +28,9 @@ namespace MathNet.Numerics.Distributions { using System; using System.Collections.Generic; + using System.Linq; + using LinearAlgebra.Double; + using LinearAlgebra.Generic; using Properties; using Statistics; @@ -80,6 +83,8 @@ namespace MathNet.Numerics.Distributions /// /// Histogram instance /// The number of trials. + /// If any of the probabilities are negative or do not sum to one. + /// If is negative. public Multinomial(Histogram h, int n) { if (h == null) @@ -209,6 +214,116 @@ namespace MathNet.Numerics.Distributions } } + /// + /// Gets the mean of the distribution. + /// + public Vector Mean + { + get + { + return _n * (DenseVector)P; + } + } + + /// + /// Gets the variance of the distribution. + /// + public Vector Variance + { + get + { + // Do not use _p, because operations below will modify _p array. Use P or _p.Clone(). + var res = (DenseVector)P; + for (var i = 0; i < res.Count; i++) + { + res[i] *= _n * (1 - res[i]); + } + + return res; + } + } + + /// + /// Gets the skewness of the distribution. + /// + public Vector Skewness + { + get + { + // Do not use _p, because operations below will modify _p array. Use P or _p.Clone(). + var res = (DenseVector)P; + for (var i = 0; i < res.Count; i++) + { + res[i] = (1.0 - (2.0 * res[i])) / Math.Sqrt(_n * (1.0 - res[i]) * res[i]); + } + + return res; + } + } + + /// + /// Computes values of the probability mass function. + /// + /// Non-negative integers x1, ..., xk + /// The probability mass at location . + /// When is null. + /// When length of is not equal to event probabilities count. + public double Probability(int[] x) + { + if (null == x) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != _p.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "x"); + } + + if (x.Sum() == _n) + { + var coef = SpecialFunctions.Multinomial(_n, x); + var num = 1.0; + for (var i = 0; i < x.Length; i++) + { + num *= Math.Pow(_p[i], x[i]); + } + + return coef * num; + } + + return 0.0; + } + + /// + /// Computes values of the log probability mass function. + /// + /// Non-negative integers x1, ..., xk + /// The log probability mass at location . + /// When is null. + /// When length of is not equal to event probabilities count. + public double ProbabilityLn(int[] x) + { + if (null == x) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != _p.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "x"); + } + + if (x.Sum() == _n) + { + var coef = Math.Log(SpecialFunctions.Multinomial(_n, x)); + var num = x.Select((t, i) => t * Math.Log(_p[i])).Sum(); + return coef + num; + } + + return 0.0; + } + /// /// Samples one multinomial distributed random variable. /// diff --git a/src/Numerics/Distributions/Multivariate/NormalGamma.cs b/src/Numerics/Distributions/Multivariate/NormalGamma.cs index c7152a2f..501d4a60 100644 --- a/src/Numerics/Distributions/Multivariate/NormalGamma.cs +++ b/src/Numerics/Distributions/Multivariate/NormalGamma.cs @@ -329,6 +329,18 @@ namespace MathNet.Numerics.Distributions } } + /// + /// Gets the variance of the distribution. + /// + /// The mean of the distribution. + public MeanPrecisionPair Variance + { + get + { + return new MeanPrecisionPair(_precisionInvScale / (_meanScale * (_precisionShape - 1)), _precisionShape / Math.Sqrt(_precisionInvScale)); + } + } + /// /// Evaluates the probability density function for a NormalGamma distribution. /// @@ -349,17 +361,17 @@ namespace MathNet.Numerics.Distributions { if (Double.IsPositiveInfinity(_precisionInvScale) && _meanScale == 0.0) { - throw new NotImplementedException(); + throw new NotSupportedException(); } if (Double.IsPositiveInfinity(_precisionInvScale)) { - throw new NotImplementedException(); + throw new NotSupportedException(); } if (_meanScale <= 0.0) { - throw new NotImplementedException(); + throw new NotSupportedException(); } // double e = -0.5 * prec * (mean - _meanLocation) * (mean - _meanLocation) - prec * _precisionInvScale; @@ -389,17 +401,17 @@ namespace MathNet.Numerics.Distributions { if (Double.IsPositiveInfinity(_precisionInvScale) && _meanScale == 0.0) { - throw new NotImplementedException(); + throw new NotSupportedException(); } if (Double.IsPositiveInfinity(_precisionInvScale)) { - throw new NotImplementedException(); + throw new NotSupportedException(); } if (_meanScale <= 0.0) { - throw new NotImplementedException(); + throw new NotSupportedException(); } // double e = -0.5 * prec * (mean - _meanLocation) * (mean - _meanLocation) - prec * _precisionInvScale; diff --git a/src/Numerics/Distributions/Multivariate/Wishart.cs b/src/Numerics/Distributions/Multivariate/Wishart.cs index 71aa39b3..4480cde3 100644 --- a/src/Numerics/Distributions/Multivariate/Wishart.cs +++ b/src/Numerics/Distributions/Multivariate/Wishart.cs @@ -214,6 +214,27 @@ namespace MathNet.Numerics.Distributions } } + /// + /// Gets the variance of the distribution. + /// + /// The variance of the distribution. + public Matrix Variance + { + get + { + var res = _s.CreateMatrix(_s.RowCount, _s.ColumnCount); + for (var i = 0; i < res.RowCount; i++) + { + for (var j = 0; j < res.ColumnCount; j++) + { + res.At(i, j, _nu * ((_s.At(i, j) * _s.At(i, j)) + (_s.At(i, i) * _s.At(j, j)))); + } + } + + return res; + } + } + /// /// Evaluates the probability density function for the Wishart distribution. /// diff --git a/src/Numerics/LinearAlgebra/IO/MatrixWriter.cs b/src/Numerics/LinearAlgebra/IO/MatrixWriter.cs index 1eec21bf..e955034b 100644 --- a/src/Numerics/LinearAlgebra/IO/MatrixWriter.cs +++ b/src/Numerics/LinearAlgebra/IO/MatrixWriter.cs @@ -36,7 +36,7 @@ namespace MathNet.Numerics.LinearAlgebra.IO /// /// Base class to write a single to a file or stream. /// - /// The data type of the Matrix. It can be either: double, float, Complex, or Complex32. + /// The data type of the matrix. public abstract class MatrixWriter where TDataType : struct, IEquatable, IFormattable { /// diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index 405a546a..2339be6d 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -524,15 +524,9 @@ LinearAlgebra\Double\IO\DelimitedReader.cs - - LinearAlgebra\Double\IO\DelimitedWriter.cs - LinearAlgebra\Double\IO\MatrixReader.cs - - LinearAlgebra\Double\IO\MatrixWriter.cs - LinearAlgebra\Double\Solvers\Iterative\BiCgStab.cs diff --git a/src/UnitTests/DistributionTests/Discrete/ZipfTests.cs b/src/UnitTests/DistributionTests/Discrete/ZipfTests.cs index caba893f..6e7d1400 100644 --- a/src/UnitTests/DistributionTests/Discrete/ZipfTests.cs +++ b/src/UnitTests/DistributionTests/Discrete/ZipfTests.cs @@ -28,7 +28,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -/*namespace MathNet.Numerics.UnitTests.DistributionTests.Discrete +namespace MathNet.Numerics.UnitTests.DistributionTests.Discrete { using System; using System.Linq; @@ -165,7 +165,7 @@ } [Test] - [ExpectedException(typeof(NotImplementedException))] + [ExpectedException(typeof(NotSupportedException))] public void ValidateMedian() { var d = new Zipf(1.0, 5); @@ -229,33 +229,39 @@ [Test] public void CanSample() { - var d = new Zipf(1.0, 5); + var d = new Zipf(0.7, 5); var s = d.Sample(); + Assert.Between(s, 0, 5); } [Test] public void CanSampleSequence() { - var d = new Zipf(1.0, 5); + var d = new Zipf(0.7, 5); var ied = d.Samples(); - var e = ied.Take(5).ToArray(); + var e = ied.Take(1000).ToArray(); + foreach (var i in e) + { + Assert.Between(i, 0, 5); + } } [Test] - [Row(0.1, 1, 1.1, 1.0)] - [Row(0.1, 20, 1.1, 0.061588204519703309)] - [Row(0.1, 50, 1.1, 0.026806743865513603)] - [Row(1.0, 1, 1.1, 1.0)] - [Row(1.0, 20, 1.1, 0.27795229652440168)] - [Row(1.0, 50, 1.1, 0.22226147170498)] - [Row(0.1, 20, 15.0, 0.061588204519703309)] - [Row(0.1, 50, 15.0, 0.026806743865513603)] - [Row(1.0, 20, 15.0, 0.27795229652440168)] - [Row(1.0, 50, 15.0, 0.22226147170498)] - public void ValidateCumulativeDistribution(double s, int n, double x, double cdf) + [Row(0.1, 1, 2)] + [Row(0.1, 20, 2)] + [Row(0.1, 50, 2)] + [Row(1.0, 1, 2)] + [Row(1.0, 20, 2)] + [Row(1.0, 50, 2)] + [Row(0.1, 20, 15)] + [Row(0.1, 50, 15)] + [Row(1.0, 20, 15)] + [Row(1.0, 50, 15)] + public void ValidateCumulativeDistribution(double s, int n, int x) { var d = new Zipf(s, n); + var cdf = SpecialFunctions.GeneralHarmonic(x, s) / SpecialFunctions.GeneralHarmonic(n, s); AssertHelpers.AlmostEqual(cdf, d.CumulativeDistribution(x), 14); } } -}*/ \ No newline at end of file +} \ No newline at end of file diff --git a/src/UnitTests/DistributionTests/Multivariate/DirichletTests.cs b/src/UnitTests/DistributionTests/Multivariate/DirichletTests.cs index bc693726..b27ca104 100644 --- a/src/UnitTests/DistributionTests/Multivariate/DirichletTests.cs +++ b/src/UnitTests/DistributionTests/Multivariate/DirichletTests.cs @@ -31,6 +31,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Multivariate { using System; + using System.Linq; using MbUnit.Framework; using Distributions; @@ -165,6 +166,36 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Multivariate } } + [Test] + [Row(new[] { 0.01, 0.03, 0.5 }, 1335.32600710379)] + [Row(new[] { 0.1, 0.2, 0.3, 0.4 }, 59.1446044600076)] + public void ValidateDensity(double[] x, double res) + { + var d = new Dirichlet(new[] { 0.1, 0.3, 0.5, 0.8 }); + AssertHelpers.AlmostEqual(res, d.Density(x), 12); + } + + [Test] + [Row(new[] { 0.01, 0.03, 0.5 })] + [Row(new[] { 0.1, 0.2, 0.3, 0.4 })] + public void ValidateDensityLn(double[] x) + { + var d = new Dirichlet(new[] { 0.1, 0.3, 0.5, 0.8 }); + AssertHelpers.AlmostEqual(d.DensityLn(x), Math.Log(d.Density(x)), 12); + } + + [Test] + [Row(new[] { 0.1, 0.3, 0.5, 0.8 })] + [Row(new[] { 0.1, 0.2, 0.3, 0.4 })] + public void ValidateEntropy(double[] x) + { + var d = new Dirichlet(x); + + var sum = x.Sum(t => (t - 1) * SpecialFunctions.DiGamma(t)); + var res = SpecialFunctions.GammaLn(x.Sum()) + (x.Sum() - x.Length) * SpecialFunctions.DiGamma(x.Sum()) - sum; + AssertHelpers.AlmostEqual(res, d.Entropy, 12); + } + [Test] public void CanSampleSymmetricDirichlet() { diff --git a/src/UnitTests/DistributionTests/Multivariate/InverseWishartTests.cs b/src/UnitTests/DistributionTests/Multivariate/InverseWishartTests.cs index 1a96d9a5..7a2929d8 100644 --- a/src/UnitTests/DistributionTests/Multivariate/InverseWishartTests.cs +++ b/src/UnitTests/DistributionTests/Multivariate/InverseWishartTests.cs @@ -193,6 +193,49 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Multivariate } } + [Test, MultipleAsserts] + [Row(1.0, 2)] + [Row(2.0, 2)] + [Row(5.0, 2)] + [Row(1.0, 5)] + [Row(2.0, 5)] + [Row(5.0, 5)] + public void ValidateMode(double nu, int order) + { + var d = new InverseWishart(nu, MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order)); + + var mode = d.Mode; + for (var i = 0; i < d.S.RowCount; i++) + { + for (var j = 0; j < d.S.ColumnCount; j++) + { + Assert.AreEqual(d.S[i, j] * (1.0 / (nu + d.S.RowCount + 1.0)), mode[i, j]); + } + } + } + + [Test, MultipleAsserts] + [Row(1.0, 2)] + [Row(2.0, 2)] + [Row(5.0, 2)] + [Row(1.0, 5)] + [Row(2.0, 5)] + [Row(5.0, 5)] + public void ValidateVariance(double nu, int order) + { + var d = new InverseWishart(nu, MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order)); + + var variance = d.Variance; + for (var i = 0; i < d.S.RowCount; i++) + { + for (var j = 0; j < d.S.ColumnCount; j++) + { + var num1 = (nu - d.S.RowCount + 1) * d.S[i, j] * d.S[i, j] + (nu - d.S.RowCount - 1) * d.S[i, i] * d.S[j, j]; + var num2 = (nu - d.S.RowCount) * (nu - d.S.RowCount - 1) * (nu - d.S.RowCount - 1) * (nu - d.S.RowCount - 3); + Assert.AreEqual(num1 / num2, variance[i, j]); + } + } + } [Test] [Row(1.0, 0.03228684517430723)] [Row(2.0, 0.018096748360719193)] diff --git a/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs b/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs index a3e13ac2..410ddcf5 100644 --- a/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs +++ b/src/UnitTests/DistributionTests/Multivariate/MultinomialTests.cs @@ -31,6 +31,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Multivariate { using System; + using System.Linq; using MbUnit.Framework; using Distributions; using Statistics; @@ -103,6 +104,65 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Multivariate Assert.AreEqual("Multinomial(Dimension = 3, Number of Trails = 4)", b.ToString()); } + [Test] + [Row(new[] { 0.3, 0.7 }, 5, new[] { 0.390360029179413, -0.390360029179413 })] + [Row(new[] { 0.1, 0.3, 0.6 }, 10, new[] { 0.843274042711568, 0.276026223736942, -0.129099444873581 })] + [Row(new[] { 0.15, 0.35, 0.3, 0.2 }, 20, new[] { 0.438357003759605, 0.140642169281549, 0.195180014589707, 0.335410196624968 })] + public void ValidateSkewness(double[] p, int n, double[] res) + { + var b = new Multinomial(p, n); + for (var i = 0; i < b.P.Length; i++) + { + AssertHelpers.AlmostEqual(res[i], b.Skewness[i], 12); + } + } + + [Test] + [Row(new[] { 0.3, 0.7 }, 5, new[] { 1.05, 1.05 })] + [Row(new[] { 0.1, 0.3, 0.6 }, 10, new[] { 0.9, 2.1, 2.4 })] + [Row(new[] { 0.15, 0.35, 0.3, 0.2 }, 20, new[] { 2.55, 4.55, 4.2, 3.2 })] + public void ValidateVariance(double[] p, int n, double[] res) + { + var b = new Multinomial(p, n); + for (var i = 0; i < b.P.Length; i++) + { + AssertHelpers.AlmostEqual(res[i], b.Variance[i], 12); + } + } + + [Test] + [Row(new[] { 0.3, 0.7 }, 5, new[] { 1.5, 3.5 })] + [Row(new[] { 0.1, 0.3, 0.6 }, 10, new[] { 1.0, 3.0, 6.0 })] + [Row(new[] { 0.15, 0.35, 0.3, 0.2 }, 20, new[] {3.0, 7.0, 6.0, 4.0 })] + public void ValidateMean(double[] p, int n, double[] res) + { + var b = new Multinomial(p, n); + for (var i = 0; i < b.P.Length; i++) + { + AssertHelpers.AlmostEqual(res[i], b.Mean[i], 12); + } + } + + [Test] + [Row(new[] { 0.3, 0.7 }, new[] { 1, 9 }, 0.121060821)] + [Row(new[] { 0.1, 0.3, 0.6 }, new[] { 1, 3, 6 }, 0.105815808)] + [Row(new[] { 0.15, 0.35, 0.3, 0.2 }, new[] { 1, 1, 1, 7 }, 0.000145152)] + public void ValidateProbability(double[] p, int[] x, double res) + { + var b = new Multinomial(p, x.Sum()); + AssertHelpers.AlmostEqual(b.Probability(x), res, 12); + } + + [Test] + [Row(new[] { 1, 2, 3, 4, 5, 6, 7, 8, 9 })] + [Row(new[] { 1, 1, 1, 2, 2, 2, 3, 3, 3 })] + [Row(new[] { 5, 6, 7, 8, 7, 6, 5, 4, 3 })] + public void ValidateProbabilityLn(int[] x) + { + var b = new Multinomial(largeP, x.Sum()); + AssertHelpers.AlmostEqual(b.ProbabilityLn(x), Math.Log(b.Probability(x)), 12); + } + [Test] public void CanSetProbability() { diff --git a/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs b/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs index 5fa03450..6e279dbf 100644 --- a/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs +++ b/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs @@ -212,7 +212,19 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Multivariate ng.RandomSource = new Random(); } - /// + [Test, MultipleAsserts] + [Row(0.0, 1.0, 1.0, 1.0)] + [Row(10.0, 1.0, 2.0, 2.0)] + public void ValidateVariance(double meanLocation, double meanScale, double precShape, double precInvScale) + { + NormalGamma ng = new NormalGamma(meanLocation, meanScale, precShape, precInvScale); + var X = precInvScale / (meanScale * (precShape - 1)); + var T = precShape / Math.Sqrt(precInvScale); + Assert.AreEqual(X, ng.Variance.Mean); + Assert.AreEqual(T, ng.Variance.Precision); + } + + /// /// Test the method which samples one variable at a time. /// [Test] diff --git a/src/UnitTests/DistributionTests/Multivariate/WishartTests.cs b/src/UnitTests/DistributionTests/Multivariate/WishartTests.cs index 8ea06f68..3424e6f3 100644 --- a/src/UnitTests/DistributionTests/Multivariate/WishartTests.cs +++ b/src/UnitTests/DistributionTests/Multivariate/WishartTests.cs @@ -214,6 +214,27 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Multivariate } } + [Test, MultipleAsserts] + [Row(1.0, 2)] + [Row(2.0, 2)] + [Row(5.0, 2)] + [Row(1.0, 5)] + [Row(2.0, 5)] + [Row(5.0, 5)] + public void ValidateVariance(double nu, int order) + { + var d = new Wishart(nu, MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order)); + + var variance = d.Variance; + for (var i = 0; i < d.S.RowCount; i++) + { + for (var j = 0; j < d.S.ColumnCount; j++) + { + Assert.AreEqual(nu * (d.S[i, j] * d.S[i, j] + d.S[i, i] * d.S[j, j]), variance[i, j]); + } + } + } + [Test] [Row(1.0, 0.014644982561926487)] [Row(2.0, 0.041042499311949421)]