From e0ba8f3aa7c70e5a9eceaed0e75eae60c800706c Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Wed, 24 Mar 2010 07:00:02 +0800 Subject: [PATCH] Added NormalGamma distribution. --- .../Distributions/Multivariate/NormalGamma.cs | 495 ++++++++++++++++++ src/Silverlight/Silverlight.csproj | 3 + .../Continuous/StudentTTests.cs | 28 +- .../Multivariate/NormalGammaTests.cs | 60 +++ 4 files changed, 574 insertions(+), 12 deletions(-) create mode 100644 src/Numerics/Distributions/Multivariate/NormalGamma.cs create mode 100644 src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs diff --git a/src/Numerics/Distributions/Multivariate/NormalGamma.cs b/src/Numerics/Distributions/Multivariate/NormalGamma.cs new file mode 100644 index 00000000..eff86d00 --- /dev/null +++ b/src/Numerics/Distributions/Multivariate/NormalGamma.cs @@ -0,0 +1,495 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.Distributions +{ + using System; + using System.Collections.Generic; + using Properties; + + /// + /// This structure represents the type over which the distribution + /// is defined. + /// + public struct MeanPrecisionPair + { + private double mMean; + private double mPrecision; + + /// + /// Constructs a new mean precision pair. + /// + /// The mean of the pair. + /// The precision of the pair. + public MeanPrecisionPair(double m, double p) + { + mMean = m; + mPrecision = p; + } + + /// + /// Gets/sets the mean of the pair. + /// + public double Mean + { + get { return mMean; } + set { mMean = value; } + } + + /// + /// Gets/sets the precision of the pair. + /// + public double Precision + { + get { return mPrecision; } + set { mPrecision = value; } + } + } + + /// + /// The distribution is the conjugate prior distribution for the + /// distribution. It specifies a prior over the mean and precision of the distribution. + /// It is parameterized by four numbers: the mean location, the mean scale, the precision shape and the + /// precision inverse scale. + /// The distribution NG(mu, tau | mloc,mscale,psscale,pinvscale) = Normal(mu | mloc, 1/(mscale*tau)) * Gamma(tau | psscale,pinvscale). + /// The following degenerate cases are special: when the precision is known, + /// the precision shape will encode the value of the precision while the precision inverse scale is positive + /// infinity. When the mean is known, the mean location will encode the value of the mean while the scale + /// will be positive infinity. A completely degenerate NormalGamma distribution with known mean and precision is possible as well. + /// + /// 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 NormalGamma + { + /// + /// The location of the mean. + /// + private double _meanLocation; + + /// + /// The scale of the mean. + /// + private double _meanScale; + + /// + /// The shape of the precision. + /// + private double _precisionShape; + + /// + /// The inverse scale of the precision. + /// + private double _precisionInvScale; + + /// + /// The distribution's random number generator. + /// + private Random _random; + + /// + /// Constructs a NormalGamma distribution. + /// + /// The location of the mean. + /// The scale of the mean. + /// The shape of the precision. + /// The inverse scale of the precision. + public NormalGamma(double meanLocation, double meanScale, double precShape, double precInvScale) + { + SetParameters(meanLocation, meanScale, precShape, precInvScale); + _random = new Random(); + } + + /// + /// Checks whether the parameters of the distribution are valid. + /// + /// The location of the mean. + /// The scale of the mean. + /// The shape of the precision. + /// The inverse scale of the precision. + /// True when the parameters are valid, false otherwise. + private static bool IsValidParameterSet(double meanLocation, double meanScale, double precShape, double precInvScale) + { + if (meanScale <= 0.0 || precShape <= 0.0 || precInvScale <= 0.0 + || Double.IsNaN(meanLocation) || Double.IsNaN(meanScale) || Double.IsNaN(precShape) + || Double.IsNaN(precInvScale)) + { + return false; + } + + return true; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// The location of the mean. + /// The scale of the mean. + /// The shape of the precision. + /// The inverse scale of the precision. + /// When the parameters don't pass the function. + private void SetParameters(double meanLocation, double meanScale, double precShape, double precInvScale) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(meanLocation, meanScale, precShape, precInvScale)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _meanLocation = meanLocation; + _meanScale = meanScale; + _precisionShape = precShape; + _precisionInvScale = precInvScale; + } + + /// + /// A string representation of the distribution. + /// + public override string ToString() + { + return "NormalGamma(Mean Location = " + _meanLocation + ", Mean Scale = " + _meanScale + + ", Precision Shape = " + _precisionShape + ", Precision Inverse Scale = " + _precisionInvScale + ")"; + } + + /// + /// Gets the location of the mean. + /// + public double MeanLocation + { + get { return _meanLocation; } + } + + /// + /// Gets the scale of the mean. + /// + public double MeanScale + { + get { return _meanScale; } + } + + /// + /// Gets the shape of the precision. + /// + public double PrecisionShape + { + get { return _precisionShape; } + } + + /// + /// Gets the inverse scale of the precision. + /// + public double PrecisionInverseScale + { + get { return _precisionInvScale; } + } + + /// + /// Returns the marginal distribution for the mean of the distribution. + /// + /// + public StudentT MeanMarginal() + { + return new StudentT(_meanLocation, _meanScale * _precisionShape / _precisionInvScale, 2.0 * _precisionShape); + } + + /// + /// Returns the marginal distribution for the precision of the distribution. + /// + /// + public Gamma PrecisionMarginal() + { + return new Gamma(_precisionShape, _precisionInvScale); + } + + /* + /// + /// Gets the mean of the distribution. + /// + /// The mean of the distribution. + public MeanPrecisionPair Mean + { + get + { + if (Double.IsPositiveInfinity(_precisionInvScale)) + { + return new MeanPrecisionPair(_meanLocation, _precisionShape); + } + else + { + return new MeanPrecisionPair(_meanLocation, _precisionShape / _precisionInvScale); + } + } + } + + /// + /// Gets or sets the random number generator. + /// + /// The random number generator used to generate a random sample. + public System.Random RandomNumberGenerator { get; set; } + + /// + /// The mode of the distribution. + /// + /// + public MeanPrecisionPair Mode + { + get + { + if (Double.IsPositiveInfinity(_precisionInvScale)) + { + return new MeanPrecisionPair(_meanLocation, _precisionShape); + } + else + { + return new MeanPrecisionPair(_meanLocation, _precisionShape / _precisionInvScale); + } + } + } + + /// + /// The median of the distribution. + /// + /// + public MeanPrecisionPair Median + { + get + { + if (Double.IsPositiveInfinity(_precisionInvScale)) + { + return new MeanPrecisionPair(_meanLocation, _precisionShape); + } + else + { + return new MeanPrecisionPair(_meanLocation, _precisionShape / _precisionInvScale); + } + } + } + + /// + /// Evaluates the probability density function for a NormalGamma distribution. + /// + public double Density(MeanPrecisionPair mp) + { + return Density(mp.Mean, mp.Precision); + } + + /// + /// Evaluates the probability density function for a NormalGamma distribution. + /// + public double Density(double mean, double prec) + { + if (Double.IsPositiveInfinity(_precisionInvScale) && _meanScale == 0.0) + { + throw new NotImplementedException(); + } + else if (Double.IsPositiveInfinity(_precisionInvScale)) + { + throw new NotImplementedException(); + } + else if (_meanScale == 0.0) + { + throw new NotImplementedException(); + } + else + { + double e = -0.5 * prec * (mean - _meanLocation) * (mean - _meanLocation) - prec * _precisionInvScale; + return System.Math.Pow(prec * _precisionInvScale, _precisionShape) * System.Math.Exp(e) + / (Math.Constants.Sqrt2Pi * System.Math.Sqrt(prec) * Math.SpecialFunctions.Gamma(_precisionShape)); + } + } + + /// + /// Evaluates the log probability density function for a NormalGamma distribution. + /// + public double DensityLn(MeanPrecisionPair mp) + { + return DensityLn(mp.Mean, mp.Precision); + } + + /// + /// Evaluates the log probability density function for a NormalGamma distribution. + /// + public double DensityLn(double mean, double prec) + { + if (Double.IsPositiveInfinity(_precisionInvScale) && _meanScale == 0.0) + { + throw new NotImplementedException(); + } + else if (Double.IsPositiveInfinity(_precisionInvScale)) + { + throw new NotImplementedException(); + } + else if (_meanScale == 0.0) + { + throw new NotImplementedException(); + } + else + { + double e = -0.5 * prec * (mean - _meanLocation) * (mean - _meanLocation) - prec * _precisionInvScale; + return (_precisionShape - 0.5) * System.Math.Log(prec) + _precisionShape * System.Math.Log(_precisionInvScale) + e + - Math.Constants.LogSqrt2Pi - Math.SpecialFunctions.GammaLn(_precisionShape); + } + } + + /// + /// Samples a NormalGamma distributed random variable. + /// + /// A random number from this distribution. + public MeanPrecisionPair Sample() + { + return NormalGamma.Sample(RandomNumberGenerator, _meanLocation, _meanScale, _precisionShape, _precisionInvScale); + } + + /// + /// Samples an array of NormalGamma distributed random variables. + /// + /// The number of variables needed. + /// An array of random numbers from this distribution. + public MeanPrecisionPair[] Sample(int size) + { + return NormalGamma.Sample(RandomNumberGenerator, size, _meanLocation, _meanScale, _precisionShape, _precisionInvScale); + } + + /// + /// Checks the parameters of a NormalGamma distribution. + /// + /// The scale of the mean. + /// The shape of the precision. + /// The inverse scale of the precision. + /// If the mean scale is negative. + /// If the inverse precision scale is negative. + /// If the precision shape is negative. + private static void CheckParameters(double meanScale, double precShape, double precInvScale) + { + if (meanScale < 0.0) + { + throw new ArgumentOutOfRangeException("meanScale", Resources.ParameterCannotBeNegative); + } + else if (precShape <= 0.0) + { + throw new ArgumentOutOfRangeException("precShape", Resources.ParameterCannotBeNegative); + } + else if (precInvScale <= 0.0) + { + throw new ArgumentOutOfRangeException("precInvScale", Resources.ParameterCannotBeNegative); + } + } + + /// + /// Samples an array of NormalGamma distributed random variables. + /// + /// The random number generator to use. + /// The location of the mean. + /// The scale of the mean. + /// The shape of the precision. + /// The inverse scale of the precision. + public static MeanPrecisionPair Sample(System.Random rnd, double meanLocation, double meanScale, double precShape, double precInvScale) + { + if (Control.CheckDistributionParameters) + { + CheckParameters(meanScale, precShape, precInvScale); + } + + MeanPrecisionPair mp = new MeanPrecisionPair(); + + // Sample the precision. + if(Double.IsPositiveInfinity(precInvScale)) + { + mp.Precision = precShape; + } + else + { + mp.Precision = Gamma.Sample(rnd, precShape, precInvScale); + } + + // Sample the mean. + if (meanScale == 0.0) + { + mp.Mean = meanLocation; + } + else + { + mp.Mean = Normal.Sample(rnd, meanLocation, System.Math.Sqrt(meanScale / mp.Precision)); + } + + return mp; + } + + /// + /// Samples an array of NormalGamma distributed random variables. + /// + /// The random number generator to use. + /// The number of variables needed. + /// The location of the mean. + /// The scale of the mean. + /// The shape of the precision. + /// The inverse scale of the precision. + public static MeanPrecisionPair[] Sample(System.Random rnd, int n, double meanLocation, double meanScale, double precShape, double precInvScale) + { + if (Control.CheckDistributionParameters) + { + CheckParameters(meanScale, precShape, precInvScale); + } + + // First sample all the precisions independently. + double[] precs = null; + if (Double.IsPositiveInfinity(precInvScale)) + { + precs = new double[n]; + for (int i = 0; i < n; i++) + { + precs[i] = precShape; + } + } + else + { + precs = Gamma.Sample(rnd, n, precShape, precInvScale); + } + + // Construct all the mean precision pairs. + MeanPrecisionPair[] arr = new MeanPrecisionPair[n]; + + // Conditionally sample all the mean. + for (int i = 0; i < n; i++) + { + arr[i].Precision = precs[i]; + if (meanScale == 0.0) + { + arr[i].Mean = meanLocation; + } + else + { + arr[i].Mean = Normal.Sample(rnd, meanLocation, System.Math.Sqrt(meanScale / precs[i])); + } + } + + return arr; + }*/ + } +} \ No newline at end of file diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index 068a9b4c..292be4ad 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -128,6 +128,9 @@ Distributions\Multivariate\Multinomial.cs + + Distributions\Multivariate\NormalGamma.cs + GlobalizationHelper.cs diff --git a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs index e6571239..c6a0ed3d 100644 --- a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs @@ -51,21 +51,25 @@ namespace MathNet.Numerics.UnitTests.DistributionTests AssertEx.AreEqual(1.0, n.DegreesOfFreedom); } - /*[Test, MultipleAsserts] - [Row(0.0, 0.0)] - [Row(0.0, 0.1)] - [Row(0.0, 1.0)] - [Row(0.0, 10.0)] - [Row(10.0, 1.0)] - [Row(-5.0, 100.0)] - [Row(0.0, Double.PositiveInfinity)] - public void CanCreateNormal(double mean, double sdev) + [Test, MultipleAsserts] + [Row(0.0, 1.0, 1.0)] + [Row(0.0, 0.1, 1.0)] + [Row(0.0, 1.0, 1.0)] + [Row(0.0, 10.0, 1.0)] + [Row(0.0, 10.0, Double.PositiveInfinity)] + [Row(10.0, 1.0, 1.0)] + [Row(-5.0, 100.0, 1.0)] + [Row(0.0, Double.PositiveInfinity, 1.0)] + public void CanCreateStudentT(double location, double scale, double dof) { - var n = new Normal(mean, sdev); - AssertEx.AreEqual(mean, n.Mean); - AssertEx.AreEqual(sdev, n.StdDev); + var n = new StudentT(location, scale, dof); + AssertEx.AreEqual(0.0, n.Location); + AssertEx.AreEqual(1.0, n.Scale); + AssertEx.AreEqual(1.0, n.DegreesOfFreedom); } + /* + [Test] [ExpectedException(typeof(ArgumentOutOfRangeException))] [Row(Double.NaN, 1.0)] diff --git a/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs b/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs new file mode 100644 index 00000000..84ffbed9 --- /dev/null +++ b/src/UnitTests/DistributionTests/Multivariate/NormalGammaTests.cs @@ -0,0 +1,60 @@ +// +// 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 NormalGammaTests + { + [Test, MultipleAsserts] + public void NormalGammaTest() + { + NormalGamma ng = new NormalGamma(10.0, 1.0, 2.0, 2.0); + + AssertEx.AreEqual(10.0, ng.MeanLocation); + AssertEx.AreEqual(1.0, ng.MeanScale); + AssertEx.AreEqual(2.0, ng.PrecisionShape); + AssertEx.AreEqual(2.0, ng.PrecisionInverseScale); + } + + [Test] + [Row(1.0, -1.3, 2.0, 2.0)] + [Row(1.0, 1.0, -1.0, 1.0)] + [Row(1.0, 1.0, 1.0, -1.0)] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void InvalidParams(double a, double b, double c, double d) + { + var nb = new NormalGamma(a, b, c, d); + } + } +} \ No newline at end of file