From 82d386f3acdf6eace46c885f002dc71fb6ea6de5 Mon Sep 17 00:00:00 2001 From: jvangael Date: Fri, 17 Jul 2009 06:25:45 +0800 Subject: [PATCH] More unit tests. Added Erfc and ErfcInv. Fixed StyleCop problems. Signed-off-by: Christoph Ruegg --- .../Continuous/NormalTests.cs | 82 ++++++- .../Managed.UnitTests.csproj | 1 + .../SpecialFunctionsTest/ErfTests.cs | 79 +++++-- src/Managed/Constants.cs | 2 +- .../Distributions/Continuous/Normal.cs | 148 ++++++++---- .../Distributions/Discrete/Bernoulli.cs | 91 +++---- .../Distributions/IContinuousDistribution.cs | 38 ++- .../Distributions/IDiscreteDistribution.cs | 14 +- src/Managed/Distributions/IDistribution.cs | 14 +- src/Managed/SiConstants.cs | 2 +- src/Managed/SiPrefixes.cs | 2 +- src/Managed/SpecialFunctions.cs | 18 +- src/Managed/SpecialFunctions/Erf.cs | 222 +++++++++++++++++- src/Native.UnitTests/Native.UnitTests.csproj | 6 + 14 files changed, 540 insertions(+), 179 deletions(-) diff --git a/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs b/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs index fd2659fc..782bda14 100644 --- a/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs +++ b/src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -26,7 +26,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.UnitTests +namespace MathNet.Numerics.UnitTests.DistributionTests { using System; using MbUnit.Framework; @@ -98,8 +98,8 @@ namespace MathNet.Numerics.UnitTests public void CanCreateNormalFromMeanAndVariance(double mean, double var) { var n = Normal.WithMeanVariance(mean, var); - AssertEx.AreEqual(mean, n.Mean); - AssertEx.AreEqual(var, n.Variance); + AssertHelpers.AlmostEqual(mean, n.Mean, 16); + AssertHelpers.AlmostEqual(var, n.Variance, 16); } [Test, MultipleAsserts] @@ -113,8 +113,8 @@ namespace MathNet.Numerics.UnitTests public void CanCreateNormalFromMeanAndPrecision(double mean, double prec) { var n = Normal.WithMeanAndPrecision(mean, prec); - AssertEx.AreEqual(mean, n.Mean); - AssertEx.AreEqual(prec, n.Precision); + AssertHelpers.AlmostEqual(mean, n.Mean, 15); + AssertHelpers.AlmostEqual(prec, n.Precision, 15); } [Test] @@ -324,6 +324,74 @@ namespace MathNet.Numerics.UnitTests } } - test samplers + [Test] + public void CanSampleStatic() + { + var d = Normal.Sample(new Random(), 0.0, 1.0); + } + + [Test] + public void CanSampleSequenceStatic() + { + var ied = Normal.Samples(new Random(), 0.0, 1.0); + var e = ied.GetEnumerator(); + e.MoveNext(); + var d = e.Current; + e.MoveNext(); + var g = e.Current; + } + + [Test] + public void CanSample() + { + var n = new Normal(); + var d = n.Sample(); + } + + [Test] + public void CanSampleSequence() + { + var n = new Normal(); + var ied = n.Samples(); + var e = ied.GetEnumerator(); + e.MoveNext(); + var d = e.Current; + e.MoveNext(); + var g = e.Current; + } + + [Test] + [Row(Double.NegativeInfinity, 0.0)] + [Row(-5.0, 0.00000028665157187919391167375233287464535385442301361187883)] + [Row(-2.0, 0.0002326290790355250363499258867279847735487493358890356)] + [Row(-0.0, 0.0062096653257761351669781045741922211278977469230927036)] + [Row(0.0, 0.0062096653257761351669781045741922211278977469230927036)] + [Row(4.0, 0.30853753872598689636229538939166226011639782444542207)] + [Row(5.0, 0.5)] + [Row(6.0, 0.69146246127401310363770461060833773988360217555457859)] + [Row(10.0, 0.9937903346742238648330218954258077788721022530769078)] + [Row(Double.PositiveInfinity, 1.0)] + public void ValidateCumulativeDistribution(double x, double f) + { + var n = Normal.WithMeanStdDev(5.0, 2.0); + AssertHelpers.AlmostEqual(f, n.CumulativeDistribution(x), 10); + } + + [Test] + [Row(Double.NegativeInfinity, 0.0)] + [Row(-5.0, 0.00000028665157187919391167375233287464535385442301361187883)] + [Row(-2.0, 0.0002326290790355250363499258867279847735487493358890356)] + [Row(-0.0, 0.0062096653257761351669781045741922211278977469230927036)] + [Row(0.0, 0.0062096653257761351669781045741922211278977469230927036)] + [Row(4.0, 0.30853753872598689636229538939166226011639782444542207)] + [Row(5.0, 0.5)] + [Row(6.0, 0.69146246127401310363770461060833773988360217555457859)] + [Row(10.0, 0.9937903346742238648330218954258077788721022530769078)] + [Row(Double.PositiveInfinity, 1.0)] + public void ValidateInverseCumulativeDistribution(double x, double f) + { + var n = Normal.WithMeanStdDev(5.0, 2.0); + AssertHelpers.AlmostEqual(x, n.InverseCumulativeDistribution(f), 10); + } } } diff --git a/src/Managed.UnitTests/Managed.UnitTests.csproj b/src/Managed.UnitTests/Managed.UnitTests.csproj index cc88c6b1..05f7aba1 100644 --- a/src/Managed.UnitTests/Managed.UnitTests.csproj +++ b/src/Managed.UnitTests/Managed.UnitTests.csproj @@ -57,6 +57,7 @@ + diff --git a/src/Managed.UnitTests/SpecialFunctionsTest/ErfTests.cs b/src/Managed.UnitTests/SpecialFunctionsTest/ErfTests.cs index f16a9073..2ae1996e 100644 --- a/src/Managed.UnitTests/SpecialFunctionsTest/ErfTests.cs +++ b/src/Managed.UnitTests/SpecialFunctionsTest/ErfTests.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -26,43 +26,70 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.UnitTests +namespace MathNet.Numerics.UnitTests.SpecialFunctionTests { - using System; - using System.IO; - using System.Collections.Generic; using MbUnit.Framework; using MathNet.Numerics; [TestFixture] public class ErfTests { - private List mLargePrecisionVals; - - [SetUp] - public void ReadLargePrecisionValues() + [Test] + [Row(0.0, 0.0)] + [Row(0.1, 0.1124629160182848984047122510143040617233925185058162)] + [Row(0.2, 0.22270258921047846617645303120925671669511570710081967)] + [Row(0.3, 0.32862675945912741618961798531820303325847175931290341)] + [Row(0.4, 0.42839235504666847645410962730772853743532927705981257)] + [Row(0.5, 0.5204998778130465376827466538919645287364515757579637)] + [Row(1.0, 0.84270079294971486934122063508260925929606699796630291)] + [Row(1.5, 0.96610514647531072706697626164594785868141047925763678)] + [Row(2.0, 0.99532226501895273416206925636725292861089179704006008)] + [Row(2.5, 0.99959304798255504106043578426002508727965132259628658)] + [Row(3.0, 0.99997790950300141455862722387041767962015229291260075)] + [Row(4.0, 0.99999998458274209971998114784032651311595142785474641)] + [Row(5.0, 0.99999999999846254020557196514981165651461662110988195)] + [Row(double.PositiveInfinity, 1.0)] + public void ErfCanMatchLargePrecision(double x, double f) { - var sr = new StreamReader(@"..\..\data\erf.txt"); - mLargePrecisionVals = new List(); + AssertHelpers.AlmostEqual(f, SpecialFunctions.Erf(x), 15); + } - while(!sr.EndOfStream) - { - var line = sr.ReadLine(); - var vals = line.Split(new char[] { '\t' }, StringSplitOptions.RemoveEmptyEntries); - mLargePrecisionVals.Add(new double[] { Double.Parse(vals[0]), Double.Parse(vals[1]) }); - } - sr.Close(); + [Test] + [Row(0.0, 1.0)] + [Row(0.1, 0.88753708398171510159528774898569593827660748149418343)] + [Row(0.2, 0.77729741078952153382354696879074328330488429289918085)] + [Row(0.3, 0.67137324054087258381038201468179696674152824068709621)] + [Row(0.4, 0.57160764495333152354589037269227146256467072294018715)] + [Row(0.5, 0.47950012218695346231725334610803547126354842424203654)] + [Row(1.0, 0.15729920705028513065877936491739074070393300203369719)] + [Row(1.5, 0.033894853524689272933023738354052141318589520742363247)] + [Row(2.0, 0.0046777349810472658379307436327470713891082029599399245)] + [Row(2.5, 0.00040695201744495893956421573997491272034867740371342016)] + [Row(3.0, 0.00002209049699858544137277612958232037984770708739924966)] + [Row(4.0, 0.000000015417257900280018852159673486884048572145253589191167)] + [Row(5.0, 0.0000000000015374597944280348501883434853833788901180503147233804)] + [Row(double.PositiveInfinity, 0.0)] + public void ErfcCanMatchLargePrecision(double x, double f) + { + AssertHelpers.AlmostEqual(f, SpecialFunctions.Erfc(x), 13); } - [Test, MultipleAsserts] - public void CanMatchLargePrecision() + [Test] + [Row(0.0, double.PositiveInfinity)] + [Row(1e-100, 15.065574702593)] // From dnA tests. + [Row(1e-30, 8.1486162231699)] // From dnA tests. + [Row(1e-20, 6.6015806223551)] // From dnA tests. + [Row(1e-10, 4.5728249585449249378479309946884581365517663258840893)] + [Row(1e-5, 3.1234132743415708640270717579666062107939039971365252)] + [Row(0.1, 1.1630871536766741628440954340547000483801487126688552)] + [Row(0.2, 0.90619380243682330953597079527631536107443494091638384)] + [Row(0.5, 0.47693627620446987338141835364313055980896974905947083)] + [Row(1.0, 0.0)] + [Row(1.5, -0.47693627620446987338141835364313055980896974905947083)] + [Row(2.0, double.NegativeInfinity)] + public void ErfcInvCanMatchLargePrecision(double x, double f) { - foreach (var xf in mLargePrecisionVals) - { - double x = xf[0]; - double f = xf[1]; - AssertEx.AreEqual(f, SpecialFunctions.Erf(x)); - } + AssertHelpers.AlmostEqual(f, SpecialFunctions.ErfcInv(x), 8); } } } diff --git a/src/Managed/Constants.cs b/src/Managed/Constants.cs index d8e62118..ab358298 100644 --- a/src/Managed/Constants.cs +++ b/src/Managed/Constants.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // diff --git a/src/Managed/Distributions/Continuous/Normal.cs b/src/Managed/Distributions/Continuous/Normal.cs index 5318ee32..6430686e 100644 --- a/src/Managed/Distributions/Continuous/Normal.cs +++ b/src/Managed/Distributions/Continuous/Normal.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -37,13 +37,18 @@ namespace MathNet.Numerics.Distributions /// public class Normal : IContinuousDistribution { - // Keeps track of the mean of the normal distribution. + /// + /// Keeps track of the mean of the normal distribution. + /// private double mMean; - // Keeps track of the standard deviation of the normal distribution. + + /// + /// Keeps track of the standard deviation of the normal distribution. + /// private double mStdDev; /// - /// Constructs a standard normal distribution. This is a normal distribution with mean 0.0 + /// Initializes a new instance of the Normal class. This is a normal distribution with mean 0.0 /// and standard deviation 1.0. The distribution will /// be initialized with the default random number generator. /// @@ -52,7 +57,7 @@ namespace MathNet.Numerics.Distributions } /// - /// Construct a normal distribution with a particular mean and standard deviation. The distribution will + /// Initializes a new instance of the Normal class with a particular mean and standard deviation. The distribution will /// be initialized with the default random number generator. /// /// The mean of the normal distribution. @@ -69,6 +74,7 @@ namespace MathNet.Numerics.Distributions /// /// The mean of the normal distribution. /// The standard deviation of the normal distribution. + /// a normal distribution. public static Normal WithMeanStdDev(double mean, double stddev) { return new Normal(mean, stddev); @@ -79,7 +85,8 @@ namespace MathNet.Numerics.Distributions /// be initialized with the default random number generator. /// /// The mean of the normal distribution. - /// The variance of the normal distribution. + /// The variance of the normal distribution. + /// a normal distribution. public static Normal WithMeanVariance(double mean, double var) { return new Normal(mean, System.Math.Sqrt(var)); @@ -90,16 +97,17 @@ namespace MathNet.Numerics.Distributions /// be initialized with the default random number generator. /// /// The mean of the normal distribution. - /// The precision of the normal distribution. + /// The precision of the normal distribution. + /// a normal distribution. public static Normal WithMeanAndPrecision(double mean, double prec) { return new Normal(mean, 1.0 / System.Math.Sqrt(prec)); } - /// /// A string representation of the distribution. /// + /// a string representation of the distribution. public override string ToString() { return "Normal(Mean = " + mMean + ", StdDev = " + mStdDev + ")"; @@ -149,23 +157,38 @@ namespace MathNet.Numerics.Distributions } /// - /// The precision of the normal distribution. + /// Gets or sets the precision of the normal distribution. /// public double Precision { - get { return 1.0 / (mStdDev * mStdDev); } - set { SetParameters(mMean, 1.0/Math.Sqrt(value)); } + get + { + return 1.0 / (mStdDev * mStdDev); + } + + set + { + double sdev = 1.0/Math.Sqrt(value); + + // Handle the case when the precision is -0. + if(Double.IsInfinity(sdev)) + { + sdev = Double.PositiveInfinity; + } + + SetParameters(mMean, sdev); + } } #region IDistribution implementation /// - /// The random number generator which is used to draw random samples. + /// Gets or sets the random number generator which is used to draw random samples. /// public Random RandomSource { get; set; } /// - /// The mean of the normal distribution. + /// Gets or sets the mean of the normal distribution. /// public double Mean { @@ -174,7 +197,7 @@ namespace MathNet.Numerics.Distributions } /// - /// The variance of the normal distribution. + /// Gets or sets the variance of the normal distribution. /// public double Variance { @@ -183,7 +206,7 @@ namespace MathNet.Numerics.Distributions } /// - /// The standard deviation of the normal distribution. + /// Gets or sets the standard deviation of the normal distribution. /// public double StdDev { @@ -192,42 +215,61 @@ namespace MathNet.Numerics.Distributions } /// - /// The entropy of the normal distribution. + /// Gets the entropy of the normal distribution. /// - public double Entropy { get { return Math.Log(mStdDev) + Constants.LogSqrt2PiE; } } + public double Entropy + { + get { return Math.Log(mStdDev) + Constants.LogSqrt2PiE; } + } /// - /// The skewness of the normal distribution. + /// Gets the skewness of the normal distribution. /// - public double Skewness { get { return 0.0; } } + public double Skewness + { + get { return 0.0; } + } #endregion #region IContinuousDistribution implementation /// - /// The mode of the normal distribution. + /// Gets the mode of the normal distribution. /// - public double Mode { get { return mMean; } } + public double Mode + { + get { return mMean; } + } /// - /// The median of the normal distribution. + /// Gets the median of the normal distribution. /// - public double Median { get { return mMean; } } + public double Median + { + get { return mMean; } + } /// - /// The minimum of the normal distribution. + /// Gets the minimum of the normal distribution. /// - public double Minimum { get { return System.Double.NegativeInfinity; } } + public double Minimum + { + get { return System.Double.NegativeInfinity; } + } /// - /// The maximum of the normal distribution. + /// Gets the maximum of the normal distribution. /// - public double Maximum { get { return System.Double.PositiveInfinity; } } + public double Maximum + { + get { return System.Double.PositiveInfinity; } + } /// /// Computes the density of the normal distribution. /// /// The location at which to compute the density. + /// the density at . public double Density(double x) { double d = (x - mMean) / mStdDev; @@ -238,26 +280,37 @@ namespace MathNet.Numerics.Distributions /// Computes the log density of the normal distribution. /// /// The location at which to compute the log density. + /// the log density at . public double DensityLn(double x) { double d = (x - mMean) / mStdDev; - return -0.5 * d * d - Math.Log(mStdDev) - Constants.LogSqrt2Pi; + return (-0.5 * d * d) - Math.Log(mStdDev) - Constants.LogSqrt2Pi; } - public double CumulativeDistribution(double x) { throw new NotImplementedException(); } + /// + /// Computes the cumulative distribution function of the normal distribution. + /// + /// The location at which to compute the cumulative density. + /// the cumulative density at . + public double CumulativeDistribution(double x) + { + return 0.5 * (1.0 + SpecialFunctions.Erf((x - mMean) / (mStdDev * System.Math.Sqrt(2.0)))); + } /// /// Generates a sample from the normal distribution using the Box-Muller algorithm. /// + /// a sample from the distribution. public double Sample() { double r2; - return mMean + mStdDev * SampleBoxMuller(RandomSource, out r2); + return mMean + (mStdDev * SampleBoxMuller(RandomSource, out r2)); } /// /// Generates a sequence of samples from the normal distribution using the Box-Muller algorithm. /// + /// a sequence of samples from the distribution. public IEnumerable Samples() { double r2; @@ -265,15 +318,20 @@ namespace MathNet.Numerics.Distributions while (true) { double r1 = SampleBoxMuller(RandomSource, out r2); - yield return mMean + mStdDev * r1; - yield return mMean + mStdDev * r2; + yield return mMean + (mStdDev * r1); + yield return mMean + (mStdDev * r2); } } #endregion + /// + /// Computes the inverse cumulative distribution function of the normal distribution. + /// + /// The location at which to compute the inverse cumulative density. + /// the inverse cumulative density at . public double InverseCumulativeDistribution(double p) { - throw new NotImplementedException(); + return mMean - (mStdDev * System.Math.Sqrt(2.0) * SpecialFunctions.ErfcInv(2.0 * p)); } /// @@ -282,10 +340,11 @@ namespace MathNet.Numerics.Distributions /// The random number generator to use. /// The mean of the normal distribution from which to generate samples. /// The standard deviation of the normal distribution from which to generate samples. + /// a sample from the distribution. public static double Sample(System.Random rng, double mean, double stddev) { double r2; - return mean + stddev * SampleBoxMuller(rng, out r2); + return mean + (stddev * SampleBoxMuller(rng, out r2)); } /// @@ -294,6 +353,7 @@ namespace MathNet.Numerics.Distributions /// The random number generator to use. /// The mean of the normal distribution from which to generate samples. /// The standard deviation of the normal distribution from which to generate samples. + /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rng, double mean, double stddev) { double r2; @@ -301,8 +361,8 @@ namespace MathNet.Numerics.Distributions while(true) { double r1 = SampleBoxMuller(rng, out r2); - yield return mean + stddev * r1; - yield return mean + stddev * r2; + yield return mean + (stddev * r1); + yield return mean + (stddev * r2); } } @@ -310,18 +370,20 @@ namespace MathNet.Numerics.Distributions /// Samples a pair of standard normal distributed random variables using the Box-Muller algorithm. /// /// The random number generator to use. - /// The second random number. + /// A second random number from the standard normal distribution computed as a side product. + /// a random number from the standard normal distribution. internal static double SampleBoxMuller(System.Random rnd, out double r2) { - double v1 = 2.0 * rnd.NextDouble() - 1.0; - double v2 = 2.0 * rnd.NextDouble() - 1.0; - double r = v1 * v1 + v2 * v2; + double v1 = (2.0 * rnd.NextDouble()) - 1.0; + double v2 = (2.0 * rnd.NextDouble()) - 1.0; + double r = (v1 * v1) + (v2 * v2); while (r >= 1.0 || r == 0.0) { - v1 = 2.0 * rnd.NextDouble() - 1.0; - v2 = 2.0 * rnd.NextDouble() - 1.0; - r = v1 * v1 + v2 * v2; + v1 = (2.0 * rnd.NextDouble()) - 1.0; + v2 = (2.0 * rnd.NextDouble()) - 1.0; + r = (v1 * v1) + (v2 * v2); } + double fac = System.Math.Sqrt(-2.0 * System.Math.Log(r) / r); r2 = v2 * fac; return v1 * fac; diff --git a/src/Managed/Distributions/Discrete/Bernoulli.cs b/src/Managed/Distributions/Discrete/Bernoulli.cs index f5e1e350..2188d008 100644 --- a/src/Managed/Distributions/Discrete/Bernoulli.cs +++ b/src/Managed/Distributions/Discrete/Bernoulli.cs @@ -1,64 +1,27 @@ -/*using System; -using System.Collections.Generic; -using Pnl.RandomSources; - -namespace Pnl.Distributions.Discrete -{ - public class Bernoulli : IDiscreteDistribution - { - public Bernoulli(double p) - { - throw new NotImplementedException(); - } - - - public override string ToString() - { - throw new NotImplementedException(); - } - - - private static void IsValidParameterSet(double p) - { - throw new NotImplementedException(); - } - - public void SetParameters(double p) - { - throw new NotImplementedException(); - } - - public double P - { - get { throw new NotImplementedException(); } - set { throw new NotImplementedException(); } - } - - #region IDistribution implementation - public RandomSource RandomNumberGenerator { get; set; } - - public double Mean { get { throw new NotImplementedException(); } } - public double Variance { get { throw new NotImplementedException(); } } - public double StdDev { get { throw new NotImplementedException(); } } - public double Entropy { get { throw new NotImplementedException(); } } - public double Skewness { get { throw new NotImplementedException(); } } - #endregion - - #region IContinuousDistribution implementation - public int Mode { get { throw new NotImplementedException(); } } - public int Median { get { throw new NotImplementedException(); } } - public int Minimum { get { throw new NotImplementedException(); } } - public int Maximum { get { throw new NotImplementedException(); } } - public double Probability(int k) { throw new NotImplementedException(); } - public double ProbabilityLn(int k) { throw new NotImplementedException(); } - public double CumulativeDistribution(int k) { throw new NotImplementedException(); } - - public int Sample() { throw new NotImplementedException(); } - public IEnumerable Samples() { throw new NotImplementedException(); } - #endregion - - public static int Sample(System.Random rng, double p) { throw new NotImplementedException(); } - public static IEnumerable Samples(System.Random rng, double p) { throw new NotImplementedException(); } - } -} -*/ \ No newline at end of file +// +// 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. +// \ No newline at end of file diff --git a/src/Managed/Distributions/IContinuousDistribution.cs b/src/Managed/Distributions/IContinuousDistribution.cs index 43bad7ee..c9a02805 100644 --- a/src/Managed/Distributions/IContinuousDistribution.cs +++ b/src/Managed/Distributions/IContinuousDistribution.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -35,14 +35,50 @@ namespace MathNet.Numerics.Distributions /// public interface IContinuousDistribution : IDistribution { + /// + /// Gets the mode of the distribution. + /// double Mode { get; } + + /// + /// Gets the median of the distribution. + /// double Median { get; } + + /// + /// Gets the smallest element in the domain of the distributions which can be represented by a double. + /// double Minimum { get; } + + /// + /// Gets the largest element in the domain of the distributions which can be represented by a double. + /// double Maximum { get; } + + /// + /// The probability density of the distribution. + /// + /// The location at which to compute the density. + /// the density at . double Density(double x); + + /// + /// The log probability density of the distribution. + /// + /// The location at which to compute the log density. + /// the log density at . double DensityLn(double x); + /// + /// Draws a random sample from the distribution. + /// + /// a sample from the distribution. double Sample(); + + /// + /// Draws a sequence of random samples from the distribution. + /// + /// a sequence of samples from the distribution. IEnumerable Samples(); } } diff --git a/src/Managed/Distributions/IDiscreteDistribution.cs b/src/Managed/Distributions/IDiscreteDistribution.cs index 76e27de5..ded5360b 100644 --- a/src/Managed/Distributions/IDiscreteDistribution.cs +++ b/src/Managed/Distributions/IDiscreteDistribution.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -36,22 +36,22 @@ namespace MathNet.Numerics.Distributions public interface IDiscreteDistribution : IDistribution { /// - /// The mode of the distribution. + /// Gets the mode of the distribution. /// int Mode { get; } /// - /// The median of the distribution. + /// Gets the median of the distribution. /// int Median { get; } /// - /// The smallest element in the domain of the distributions which can be represented by an integer. + /// Gets the smallest element in the domain of the distributions which can be represented by an integer. /// int Minimum { get; } /// - /// The largest element in the domain of the distributions which can be represented by an integer. + /// Gets the largest element in the domain of the distributions which can be represented by an integer. /// int Maximum { get; } @@ -59,22 +59,26 @@ namespace MathNet.Numerics.Distributions /// Computes values of the probability mass function. /// /// The location in the domain where we want to evaluate the probability mass function. + /// the probability mass at location . double Probability(int k); /// /// Computes values of the log probability mass function. /// /// The location in the domain where we want to evaluate the log probability mass function. + /// the log probability mass at location . double ProbabilityLn(int k); /// /// Draws a random sample from the distribution. /// + /// a sample from the distribution. int Sample(); /// /// Draws a sequence of random samples from the distribution. /// + /// a sequence of samples from the distribution. IEnumerable Samples(); } } diff --git a/src/Managed/Distributions/IDistribution.cs b/src/Managed/Distributions/IDistribution.cs index 86392e48..a3898b1e 100644 --- a/src/Managed/Distributions/IDistribution.cs +++ b/src/Managed/Distributions/IDistribution.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -41,33 +41,35 @@ namespace MathNet.Numerics.Distributions Random RandomSource { get; set; } /// - /// The mean of the distribution. + /// Gets the mean of the distribution. /// double Mean { get; } /// - /// The variance of the distribution. + /// Gets the variance of the distribution. /// double Variance { get; } /// - /// The standard deviation of the distribution. + /// Gets the standard deviation of the distribution. /// double StdDev { get; } /// - /// The entropy of the distribution. + /// Gets the entropy of the distribution. /// double Entropy { get; } /// - /// The skewness of the distribution. + /// Gets the skewness of the distribution. /// double Skewness { get; } /// /// Computes the cumulative distribution function (cdf) for this probability distribution. /// + /// The location at which to compute the cumulative distribution function. + /// the cumulative distribution at location . double CumulativeDistribution(double x); } } diff --git a/src/Managed/SiConstants.cs b/src/Managed/SiConstants.cs index 0582c1df..e42086f1 100644 --- a/src/Managed/SiConstants.cs +++ b/src/Managed/SiConstants.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // diff --git a/src/Managed/SiPrefixes.cs b/src/Managed/SiPrefixes.cs index aefbbd36..6f4d4e8d 100644 --- a/src/Managed/SiPrefixes.cs +++ b/src/Managed/SiPrefixes.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // diff --git a/src/Managed/SpecialFunctions.cs b/src/Managed/SpecialFunctions.cs index 6b0ada19..931c82fd 100644 --- a/src/Managed/SpecialFunctions.cs +++ b/src/Managed/SpecialFunctions.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -34,21 +34,5 @@ namespace MathNet.Numerics /// public static partial class SpecialFunctions { - /// - /// A helper function to evaluate polynomials fast. - /// - /// The coefficients of the polynomial. - /// The location where to evaluate the polynomial at. - private static double evaluate_polynomial(double[] poly, double z) - { - int count = poly.Length; - double sum = poly[count - 1]; - for (int i = count - 2; i >= 0; --i) - { - sum *= z; - sum += poly[i]; - } - return sum; - } } } diff --git a/src/Managed/SpecialFunctions/Erf.cs b/src/Managed/SpecialFunctions/Erf.cs index feab0b8e..83e42627 100644 --- a/src/Managed/SpecialFunctions/Erf.cs +++ b/src/Managed/SpecialFunctions/Erf.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -30,11 +30,14 @@ namespace MathNet.Numerics { using System; + /// + /// This partial implementation of the SpecialFunctions class contains all methods related to the error function. + /// public partial class SpecialFunctions { /// Calculates the error function. /// The value to evaluate. - /// The error function evaluated at given value. + /// the error function evaluated at given value. /// /// /// returns 1 if x == Double.PositiveInfinity. @@ -47,14 +50,17 @@ namespace MathNet.Numerics { return 0; } + if (Double.IsPositiveInfinity(x)) { return 1; } + if (Double.IsNegativeInfinity(x)) { return -1; } + if (Double.IsNaN(x) || Double.IsNaN(x)) { return Double.NaN; @@ -65,7 +71,7 @@ namespace MathNet.Numerics /// Calculates the complementary error function. /// The value to evaluate. - /// The complementary error function evaluated at given value. + /// the complementary error function evaluated at given value. /// /// /// returns 0 if x == Double.PositiveInfinity. @@ -78,14 +84,17 @@ namespace MathNet.Numerics { return 1; } + if (Double.IsPositiveInfinity(x)) { return 0; } + if (Double.IsNegativeInfinity(x)) { return 2; } + if (Double.IsNaN(x) || Double.IsNaN(x)) { return Double.NaN; @@ -94,7 +103,12 @@ namespace MathNet.Numerics return ErfImp(x,true); } - + /// + /// Implementation of the error function. + /// + /// Where to evaluate the error function. + /// Whether to compute 1 - the error function. + /// the error function. private static double ErfImp(double z, bool invert) { if (z < 0) @@ -103,10 +117,12 @@ namespace MathNet.Numerics { return -ErfImp(-z, invert); } + if (z < -0.5) { return 2 - ErfImp((-z), invert); } + return 1 + ErfImp(-z, false); } @@ -124,7 +140,7 @@ namespace MathNet.Numerics // if (z < 1e-10) { - result = z * 1.125 + z * 0.003379167095512573896158903121545171688; + result = (z * 1.125) + (z * 0.003379167095512573896158903121545171688); } else { @@ -132,7 +148,7 @@ namespace MathNet.Numerics double[] n = new double[] { 0.00337916709551257388990745, -0.00073695653048167948530905, -0.374732337392919607868241, 0.0817442448733587196071743, -0.0421089319936548595203468, 0.0070165709512095756344528, -0.00495091255982435110337458, 0.000871646599037922480317225 }; double[] d = new double[] { 1, -0.218088218087924645390535, 0.412542972725442099083918, -0.0841891147873106755410271, 0.0655338856400241519690695, -0.0120019604454941768171266, 0.00408165558926174048329689, -0.000615900721557769691924509 }; - result = z * 1.125 + z * evaluate_polynomial(n, z) / evaluate_polynomial(d, z); + result = (z * 1.125) + (z * evaluate_polynomial(n, z) / evaluate_polynomial(d, z)); } } else if ((z < 110) || ((z < 110) && invert)) @@ -246,8 +262,9 @@ namespace MathNet.Numerics r = evaluate_polynomial(n, z - 85) / evaluate_polynomial(d, z - 85); b = 0.5641584396F; } + double g = System.Math.Exp(-z * z) / z; - result = g * b + g * r; + result = (g * b) + (g * r); } else { @@ -265,5 +282,196 @@ namespace MathNet.Numerics return result; } + + ///Calculates the complementary inverse error function evaluated at z. + /// The complementary inverse error function evaluated at given value. + /// We have tested this implementation against the arbitrary precision mpmath library + /// and found cases where we can only guarantee 9 significant figures correct. + /// + /// returns Double.PositiveInfinity if z <= 0.0. + /// returns Double.NegativeInfinity if z >= 2.0. + /// + /// + ///Calculates the complementary inverse error function evaluated at z. + ///value to evaluate. + ///the complementary inverse error function evaluated at Z. + public static double ErfcInv(double z) + { + if (z <= 0.0) + { + return double.PositiveInfinity; + } + + if (z >= 2.0) + { + return double.NegativeInfinity; + } + + double p, q, s; + if (z > 1) + { + q = 2 - z; + p = 1 - q; + s = -1; + } + else + { + p = 1 - z; + q = z; + s = 1; + } + + return ErfInvImpl(p, q, s); + } + + /// + /// The implementation of the inverse error function. + /// + /// First intermediate parameter. + /// Second intermediate parameter. + /// Third intermediate parameter. + /// the inverse error function. + private static double ErfInvImpl(double p, double q, double s) + { + double result; + + if (p <= 0.5) + { + // + // Evaluate inverse erf using the rational approximation: + // + // x = p(p+10)(Y+R(p)) + // + // Where Y is a constant, and R(p) is optimized for a low + // absolute error compared to |Y|. + // + // double: Max error found: 2.001849e-18 + // long double: Max error found: 1.017064e-20 + // Maximum Deviation Found (actual error term at infinite precision) 8.030e-21 + // + float Y = 0.0891314744949340820313f; + double[] P = new double[] { -0.000508781949658280665617, -0.00836874819741736770379, 0.0334806625409744615033, -0.0126926147662974029034, -0.0365637971411762664006, 0.0219878681111168899165, 0.00822687874676915743155, -0.00538772965071242932965 }; + double[] Q = new double[] { 1, -0.970005043303290640362, -1.56574558234175846809, 1.56221558398423026363, 0.662328840472002992063, -0.71228902341542847553, -0.0527396382340099713954, 0.0795283687341571680018, -0.00233393759374190016776, 0.000886216390456424707504 }; + double g = p * (p + 10); + double r = evaluate_polynomial(P, p) / evaluate_polynomial(Q, p); + result = (g * Y) + (g * r); + } + else if (q >= 0.25) + { + // + // Rational approximation for 0.5 > q >= 0.25 + // + // x = sqrt(-2*log(q)) / (Y + R(q)) + // + // Where Y is a constant, and R(q) is optimized for a low + // absolute error compared to Y. + // + // double : Max error found: 7.403372e-17 + // long double : Max error found: 6.084616e-20 + // Maximum Deviation Found (error term) 4.811e-20 + // + float Y = 2.249481201171875f; + double[] P = new double[] { -0.202433508355938759655, 0.105264680699391713268, 8.37050328343119927838, 17.6447298408374015486, -18.8510648058714251895, -44.6382324441786960818, 17.445385985570866523, 21.1294655448340526258, -3.67192254707729348546 }; + double[] Q = new double[] { 1, 6.24264124854247537712, 3.9713437953343869095, -28.6608180499800029974, -20.1432634680485188801, 48.5609213108739935468, 10.8268667355460159008, -22.6436933413139721736, 1.72114765761200282724 }; + double g = System.Math.Sqrt(-2 * System.Math.Log(q)); + double xs = q - 0.25; + double r = evaluate_polynomial(P, xs) / evaluate_polynomial(Q, xs); + result = g / (Y + r); + } + else + { + // + // For q < 0.25 we have a series of rational approximations all + // of the general form: + // + // let: x = sqrt(-log(q)) + // + // Then the result is given by: + // + // x(Y+R(x-B)) + // + // where Y is a constant, B is the lowest value of x for which + // the approximation is valid, and R(x-B) is optimized for a low + // absolute error compared to Y. + // + // Note that almost all code will really go through the first + // or maybe second approximation. After than we're dealing with very + // small input values indeed: 80 and 128 bit long double's go all the + // way down to ~ 1e-5000 so the "tail" is rather long... + // + double x = System.Math.Sqrt(-System.Math.Log(q)); + if (x < 3) + { + // Max error found: 1.089051e-20 + float Y = 0.807220458984375f; + double[] P = new double[] { -0.131102781679951906451, -0.163794047193317060787, 0.117030156341995252019, 0.387079738972604337464, 0.337785538912035898924, 0.142869534408157156766, 0.0290157910005329060432, 0.00214558995388805277169, -0.679465575181126350155e-6, 0.285225331782217055858e-7, -0.681149956853776992068e-9 }; + double[] Q = new double[] { 1, 3.46625407242567245975, 5.38168345707006855425, 4.77846592945843778382, 2.59301921623620271374, 0.848854343457902036425, 0.152264338295331783612, 0.01105924229346489121 }; + double xs = x - 1.125; + double R = evaluate_polynomial(P, xs) / evaluate_polynomial(Q, xs); + result = (Y * x) + (R * x); + } + else if (x < 6) + { + // Max error found: 8.389174e-21 + float Y = 0.93995571136474609375f; + double[] P = new double[] { -0.0350353787183177984712, -0.00222426529213447927281, 0.0185573306514231072324, 0.00950804701325919603619, 0.00187123492819559223345, 0.000157544617424960554631, 0.460469890584317994083e-5, -0.230404776911882601748e-9, 0.266339227425782031962e-11 }; + double[] Q = new double[] { 1, 1.3653349817554063097, 0.762059164553623404043, 0.220091105764131249824, 0.0341589143670947727934, 0.00263861676657015992959, 0.764675292302794483503e-4 }; + double xs = x - 3; + double R = evaluate_polynomial(P, xs) / evaluate_polynomial(Q, xs); + result = (Y * x) + (R * x); + } + else if (x < 18) + { + // Max error found: 1.481312e-19 + float Y = 0.98362827301025390625f; + double[] P = new double[] { -0.0167431005076633737133, -0.00112951438745580278863, 0.00105628862152492910091, 0.000209386317487588078668, 0.149624783758342370182e-4, 0.449696789927706453732e-6, 0.462596163522878599135e-8, -0.281128735628831791805e-13, 0.99055709973310326855e-16 }; + double[] Q = new double[] { 1, 0.591429344886417493481, 0.138151865749083321638, 0.0160746087093676504695, 0.000964011807005165528527, 0.275335474764726041141e-4, 0.282243172016108031869e-6 }; + double xs = x - 6; + double R = evaluate_polynomial(P, xs) / evaluate_polynomial(Q, xs); + result = (Y * x) + (R * x); + } + else if (x < 44) + { + // Max error found: 5.697761e-20 + float Y = 0.99714565277099609375f; + double[] P = new double[] { -0.0024978212791898131227, -0.779190719229053954292e-5, 0.254723037413027451751e-4, 0.162397777342510920873e-5, 0.396341011304801168516e-7, 0.411632831190944208473e-9, 0.145596286718675035587e-11, -0.116765012397184275695e-17 }; + double[] Q = new double[] { 1, 0.207123112214422517181, 0.0169410838120975906478, 0.000690538265622684595676, 0.145007359818232637924e-4, 0.144437756628144157666e-6, 0.509761276599778486139e-9 }; + double xs = x - 18; + double R = evaluate_polynomial(P, xs) / evaluate_polynomial(Q, xs); + result = (Y * x) + (R * x); + } + else + { + // Max error found: 1.279746e-20 + float Y = 0.99941349029541015625f; + double[] P = new double[] { -0.000539042911019078575891, -0.28398759004727721098e-6, 0.899465114892291446442e-6, 0.229345859265920864296e-7, 0.225561444863500149219e-9, 0.947846627503022684216e-12, 0.135880130108924861008e-14, -0.348890393399948882918e-21 }; + double[] Q = new double[] { 1, 0.0845746234001899436914, 0.00282092984726264681981, 0.468292921940894236786e-4, 0.399968812193862100054e-6, 0.161809290887904476097e-8, 0.231558608310259605225e-11 }; + double xs = x - 44; + double R = evaluate_polynomial(P, xs) / evaluate_polynomial(Q, xs); + result = (Y * x) + (R * x); + } + } + + return s * result; + } + + /// + /// A helper function to evaluate polynomials fast. + /// + /// The coefficients of the polynomial. + /// The location where to evaluate the polynomial at. + /// the evaluation of the polynomial. + private static double evaluate_polynomial(double[] poly, double z) + { + int count = poly.Length; + double sum = poly[count - 1]; + for (int i = count - 2; i >= 0; --i) + { + sum *= z; + sum += poly[i]; + } + + return sum; + } } } diff --git a/src/Native.UnitTests/Native.UnitTests.csproj b/src/Native.UnitTests/Native.UnitTests.csproj index c7a2372a..3b5a099a 100644 --- a/src/Native.UnitTests/Native.UnitTests.csproj +++ b/src/Native.UnitTests/Native.UnitTests.csproj @@ -57,6 +57,9 @@ + + AssertHelpers.cs + CombinatoricsTests\CombinatoricsCountingTest.cs @@ -66,6 +69,9 @@ DistributionTests\Continuous\NormalTests.cs + + PrecisionTest.cs + SpecialFunctionsTest\ErfTests.cs