merged Andriy's distribution additions and corrections

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);
             }
         }
 
-        /// <summary>
+        /// <summary> 
         /// Computes the cumulative distribution function of the distribution.
         /// </summary>
-        /// <param name="x">The location at which to compute the cumulative density.</param>
+        /// <param name="x">The integer location at which to compute the cumulative density.</param>
         /// <returns>the cumulative density at <paramref name="x"/>.</returns>
         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
         /// <returns>a random number from the Zipf distribution.</returns>
         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;
         }
     }
 }

src/Numerics/Distributions/Multivariate/Dirichlet.cs

@@ -215,6 +215,80 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
+        /// <summary>
+        /// Gets the entropy of the distribution.
+        /// </summary>
+        public double Entropy
+        {
+            get
+            {
+                var num = _alpha.Sum(t => (t - 1) * SpecialFunctions.DiGamma(t));
+                return SpecialFunctions.GammaLn(AlphaSum) + ((AlphaSum - Dimension) * SpecialFunctions.DiGamma(AlphaSum)) - num;
+            }
+        }
+
+        /// <summary>
+        /// Computes the density of the distribution.
+        /// </summary>
+        /// <param name="x">The locations at which to compute the density.</param>
+        /// <returns>the density at <paramref name="x"/>.</returns>
+        /// <remarks>The Dirichlet distribution requires that the sum of the components of x equals 1. 
+        /// You can also leave out the last <paramref name="x"/> component, and it will be computed from the others. </remarks>
+        public double Density(double[] x)
+        {
+            return Math.Exp(DensityLn(x));
+        }
+
+        /// <summary>
+        /// Computes the log density of the distribution.
+        /// </summary>
+        /// <param name="x">The locations at which to compute the density.</param>
+        /// <returns>the density at <paramref name="x"/>.</returns>
+        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;
+        }
+
         /// <summary>
         /// Gets or sets the random number generator which is used to draw random samples.
         /// </summary>

src/Numerics/Distributions/Multivariate/InverseWishart.cs

@@ -200,6 +200,43 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
+        /// <summary>
+        /// Gets the mode of the distribution.
+        /// </summary>
+        /// <value>The mode of the distribution.</value>
+        /// <remarks>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.</remarks>
+        public Matrix<double> Mode
+        {
+            get
+            {
+                return _s * (1.0 / (_nu + _s.RowCount + 1.0));
+            }
+        }
+
+        /// <summary>
+        /// Gets the variance of the distribution.
+        /// </summary>
+        /// <value>The variance  of the distribution.</value>
+        /// <remarks>Kanti V. Mardia, J. T. Kent and J. M. Bibby (1979). Multivariate Analysis.</remarks>
+        public Matrix<double> 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;
+            }
+        }
+
         /// <summary>
         /// Evaluates the probability density function for the inverse Wishart distribution.
         /// </summary>

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
         /// </summary>
         /// <param name="h">Histogram instance</param>
         /// <param name="n">The number of trials.</param>
+                /// <exception cref="ArgumentOutOfRangeException">If any of the probabilities are negative or do not sum to one.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="n"/> is negative.</exception>
         public Multinomial(Histogram h, int n)
         {
             if (h == null)
@@ -209,6 +214,116 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
+        /// <summary>
+        /// Gets the mean of the distribution.
+        /// </summary>
+        public Vector<double> Mean
+        {
+            get
+            {
+                return _n * (DenseVector)P;
+            }
+        }
+
+        /// <summary>
+        /// Gets the variance of the distribution.
+        /// </summary>
+        public Vector<double> 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;
+            }
+        }
+
+        /// <summary>
+        /// Gets the skewness of the distribution.
+        /// </summary>
+        public Vector<double> 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;
+            }
+        }
+
+        /// <summary>
+        /// Computes values of the probability mass function.
+        /// </summary>
+        /// <param name="x">Non-negative integers x1, ..., xk</param>
+        /// <returns>The probability mass at location <paramref name="x"/>.</returns>
+        /// <exception cref="ArgumentNullException">When <paramref name="x"/> is null.</exception>
+        /// <exception cref="ArgumentException">When length of <paramref name="x"/> is not equal to event probabilities count.</exception>
+        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;
+        }
+
+        /// <summary>
+        /// Computes values of the log probability mass function.
+        /// </summary>
+        /// <param name="x">Non-negative integers x1, ..., xk</param>
+        /// <returns>The log probability mass at location <paramref name="x"/>.</returns>
+        /// <exception cref="ArgumentNullException">When <paramref name="x"/> is null.</exception>
+        /// <exception cref="ArgumentException">When length of <paramref name="x"/> is not equal to event probabilities count.</exception>
+        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;
+        }
+
         /// <summary>
         /// Samples one multinomial distributed random variable.
         /// </summary>

src/Numerics/Distributions/Multivariate/NormalGamma.cs

@@ -329,6 +329,18 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
+        /// <summary>
+        /// Gets the variance of the distribution.
+        /// </summary>
+        /// <value>The mean of the distribution.</value>
+        public MeanPrecisionPair Variance
+        {
+            get
+            {
+                return new MeanPrecisionPair(_precisionInvScale / (_meanScale * (_precisionShape - 1)), _precisionShape / Math.Sqrt(_precisionInvScale));
+            }
+        }
+
         /// <summary>
         /// Evaluates the probability density function for a NormalGamma distribution.
         /// </summary>
@@ -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;

src/Numerics/Distributions/Multivariate/Wishart.cs

@@ -214,6 +214,27 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
+        /// <summary>
+        /// Gets the variance of the distribution.
+        /// </summary>
+        /// <value>The variance  of the distribution.</value>
+        public Matrix<double> 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;
+            }
+        }
+
         /// <summary>
         /// Evaluates the probability density function for the Wishart distribution.
         /// </summary>

src/Numerics/LinearAlgebra/IO/MatrixWriter.cs

@@ -36,7 +36,7 @@ namespace MathNet.Numerics.LinearAlgebra.IO
     /// <summary>
     /// Base class to write a single <see cref="Matrix{DataType}"/> to a file or stream.
     /// </summary>
-    /// <typeparam name="TDataType">The data type of the Matrix. It can be either: double, float, Complex, or Complex32.</typeparam>
+    /// <typeparam name="TDataType">The data type of the matrix.</typeparam>
     public abstract class MatrixWriter<TDataType> where TDataType : struct, IEquatable<TDataType>, IFormattable
     {
         /// <summary>

src/Silverlight/Silverlight.csproj

@@ -524,15 +524,9 @@
     <Compile Include="..\Numerics\LinearAlgebra\Double\IO\DelimitedReader.cs">
       <Link>LinearAlgebra\Double\IO\DelimitedReader.cs</Link>
     </Compile>
-    <Compile Include="..\Numerics\LinearAlgebra\Double\IO\DelimitedWriter.cs">
-      <Link>LinearAlgebra\Double\IO\DelimitedWriter.cs</Link>
-    </Compile>
     <Compile Include="..\Numerics\LinearAlgebra\Double\IO\MatrixReader.cs">
       <Link>LinearAlgebra\Double\IO\MatrixReader.cs</Link>
     </Compile>
-    <Compile Include="..\Numerics\LinearAlgebra\Double\IO\MatrixWriter.cs">
-      <Link>LinearAlgebra\Double\IO\MatrixWriter.cs</Link>
-    </Compile>
     <Compile Include="..\Numerics\LinearAlgebra\Double\Solvers\Iterative\BiCgStab.cs">
       <Link>LinearAlgebra\Double\Solvers\Iterative\BiCgStab.cs</Link>
     </Compile>

src/UnitTests/DistributionTests/Discrete/ZipfTests.cs

@@ -28,7 +28,7 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-/*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);
         }
     }
-}*/
+}

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()
         {

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)]

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<string>("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()
         {

src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs

@@ -212,7 +212,19 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Multivariate
             ng.RandomSource = new Random();
         }
 
-        /// <summary>
+        [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);
+        }
+
+	    /// <summary>
         /// Test the method which samples one variable at a time.
         /// </summary>
         [Test]

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)]