tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

More unit tests. 

Added Erfc and ErfcInv.
Fixed StyleCop problems.

Signed-off-by: Christoph Ruegg <git@cdrnet.ch>
pull/2/head
jvangael 17 years ago
committed by Christoph Ruegg
parent
7b45d69fab
commit
82d386f3ac
14 changed files with 540 additions and 179 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 82
      src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs
  2. 1
      src/Managed.UnitTests/Managed.UnitTests.csproj
  3. 79
      src/Managed.UnitTests/SpecialFunctionsTest/ErfTests.cs
  4. 2
      src/Managed/Constants.cs
  5. 148
      src/Managed/Distributions/Continuous/Normal.cs
  6. 91
      src/Managed/Distributions/Discrete/Bernoulli.cs
  7. 38
      src/Managed/Distributions/IContinuousDistribution.cs
  8. 14
      src/Managed/Distributions/IDiscreteDistribution.cs
  9. 14
      src/Managed/Distributions/IDistribution.cs
  10. 2
      src/Managed/SiConstants.cs
  11. 2
      src/Managed/SiPrefixes.cs
  12. 18
      src/Managed/SpecialFunctions.cs
  13. 222
      src/Managed/SpecialFunctions/Erf.cs
  14. 6
      src/Native.UnitTests/Native.UnitTests.csproj

82
src/Managed.UnitTests/DistributionTests/Continuous/NormalTests.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="NormalTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
@ -26,7 +26,7 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
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<double>(mean, n.Mean);
AssertEx.AreEqual<double>(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<double>(mean, n.Mean);
AssertEx.AreEqual<double>(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);
}
}
}

1
src/Managed.UnitTests/Managed.UnitTests.csproj
View File

@ -57,6 +57,7 @@
<Reference Include="System.Xml" />
</ItemGroup>
<ItemGroup>
<Compile Include="AssertHelpers.cs" />
<Compile Include="CombinatoricsTests\CombinatoricsCountingTest.cs" />
<Compile Include="ComplexTest.cs" />
<Compile Include="DistributionTests\Continuous\NormalTests.cs" />

79
src/Managed.UnitTests/SpecialFunctionsTest/ErfTests.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="ErfTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
@ -26,43 +26,70 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
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<double []> 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<double[]>();
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<double>(f, SpecialFunctions.Erf(x));
}
AssertHelpers.AlmostEqual(f, SpecialFunctions.ErfcInv(x), 8);
}
}
}

2
src/Managed/Constants.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="Constants.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//

148
src/Managed/Distributions/Continuous/Normal.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="Normal.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
@ -37,13 +37,18 @@ namespace MathNet.Numerics.Distributions
/// </summary>
public class Normal : IContinuousDistribution
{
// Keeps track of the mean of the normal distribution.
/// <summary>
/// Keeps track of the mean of the normal distribution.
/// </summary>
private double mMean;
// Keeps track of the standard deviation of the normal distribution.
/// <summary>
/// Keeps track of the standard deviation of the normal distribution.
/// </summary>
private double mStdDev;
/// <summary>
/// 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 <seealso cref="System.Random"/> random number generator.
/// </summary>
@ -52,7 +57,7 @@ namespace MathNet.Numerics.Distributions
}
/// <summary>
/// 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 <seealso cref="System.Random"/> random number generator.
/// </summary>
/// <param name="mean">The mean of the normal distribution.</param>
@ -69,6 +74,7 @@ namespace MathNet.Numerics.Distributions
/// </summary>
/// <param name="mean">The mean of the normal distribution.</param>
/// <param name="stddev">The standard deviation of the normal distribution.</param>
/// <returns>a normal distribution.</returns>
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 <seealso cref="System.Random"/> random number generator.
/// </summary>
/// <param name="mean">The mean of the normal distribution.</param>
/// <param name="stddev">The variance of the normal distribution.</param>
/// <param name="var">The variance of the normal distribution.</param>
/// <returns>a normal distribution.</returns>
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 <seealso cref="System.Random"/> random number generator.
/// </summary>
/// <param name="mean">The mean of the normal distribution.</param>
/// <param name="stddev">The precision of the normal distribution.</param>
/// <param name="prec">The precision of the normal distribution.</param>
/// <returns>a normal distribution.</returns>
public static Normal WithMeanAndPrecision(double mean, double prec)
{
return new Normal(mean, 1.0 / System.Math.Sqrt(prec));
}
/// <summary>
/// A string representation of the distribution.
/// </summary>
/// <returns>a string representation of the distribution.</returns>
public override string ToString()
{
return "Normal(Mean = " + mMean + ", StdDev = " + mStdDev + ")";
@ -149,23 +157,38 @@ namespace MathNet.Numerics.Distributions
}
/// <summary>
/// The precision of the normal distribution.
/// Gets or sets the precision of the normal distribution.
/// </summary>
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
/// <summary>
/// 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.
/// </summary>
public Random RandomSource { get; set; }
/// <summary>
/// The mean of the normal distribution.
/// Gets or sets the mean of the normal distribution.
/// </summary>
public double Mean
{
@ -174,7 +197,7 @@ namespace MathNet.Numerics.Distributions
}
/// <summary>
/// The variance of the normal distribution.
/// Gets or sets the variance of the normal distribution.
/// </summary>
public double Variance
{
@ -183,7 +206,7 @@ namespace MathNet.Numerics.Distributions
}
/// <summary>
/// The standard deviation of the normal distribution.
/// Gets or sets the standard deviation of the normal distribution.
/// </summary>
public double StdDev
{
@ -192,42 +215,61 @@ namespace MathNet.Numerics.Distributions
}
/// <summary>
/// The entropy of the normal distribution.
/// Gets the entropy of the normal distribution.
/// </summary>
public double Entropy { get { return Math.Log(mStdDev) + Constants.LogSqrt2PiE; } }
public double Entropy
{
get { return Math.Log(mStdDev) + Constants.LogSqrt2PiE; }
}
/// <summary>
/// The skewness of the normal distribution.
/// Gets the skewness of the normal distribution.
/// </summary>
public double Skewness { get { return 0.0; } }
public double Skewness
{
get { return 0.0; }
}
#endregion
#region IContinuousDistribution implementation
/// <summary>
/// The mode of the normal distribution.
/// Gets the mode of the normal distribution.
/// </summary>
public double Mode { get { return mMean; } }
public double Mode
{
get { return mMean; }
}
/// <summary>
/// The median of the normal distribution.
/// Gets the median of the normal distribution.
/// </summary>
public double Median { get { return mMean; } }
public double Median
{
get { return mMean; }
}
/// <summary>
/// The minimum of the normal distribution.
/// Gets the minimum of the normal distribution.
/// </summary>
public double Minimum { get { return System.Double.NegativeInfinity; } }
public double Minimum
{
get { return System.Double.NegativeInfinity; }
}
/// <summary>
/// The maximum of the normal distribution.
/// Gets the maximum of the normal distribution.
/// </summary>
public double Maximum { get { return System.Double.PositiveInfinity; } }
public double Maximum
{
get { return System.Double.PositiveInfinity; }
}
/// <summary>
/// Computes the density of the normal distribution.
/// </summary>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
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.
/// </summary>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
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(); }
/// <summary>
/// Computes the cumulative distribution function of the normal distribution.
/// </summary>
/// <param name="x">The location at which to compute the cumulative density.</param>
/// <returns>the cumulative density at <paramref name="x"/>.</returns>
public double CumulativeDistribution(double x)
{
return 0.5 * (1.0 + SpecialFunctions.Erf((x - mMean) / (mStdDev * System.Math.Sqrt(2.0))));
}
/// <summary>
/// Generates a sample from the normal distribution using the <i>Box-Muller</i> algorithm.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public double Sample()
{
double r2;
return mMean + mStdDev * SampleBoxMuller(RandomSource, out r2);
return mMean + (mStdDev * SampleBoxMuller(RandomSource, out r2));
}
/// <summary>
/// Generates a sequence of samples from the normal distribution using the <i>Box-Muller</i> algorithm.
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> 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
/// <summary>
/// Computes the inverse cumulative distribution function of the normal distribution.
/// </summary>
/// <param name="p">The location at which to compute the inverse cumulative density.</param>
/// <returns>the inverse cumulative density at <paramref name="x"/>.</returns>
public double InverseCumulativeDistribution(double p)
{
throw new NotImplementedException();
return mMean - (mStdDev * System.Math.Sqrt(2.0) * SpecialFunctions.ErfcInv(2.0 * p));
}
/// <summary>
@ -282,10 +340,11 @@ namespace MathNet.Numerics.Distributions
/// <param name="rng">The random number generator to use.</param>
/// <param name="mean">The mean of the normal distribution from which to generate samples.</param>
/// <param name="stddev">The standard deviation of the normal distribution from which to generate samples.</param>
/// <returns>a sample from the distribution.</returns>
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));
}
/// <summary>
@ -294,6 +353,7 @@ namespace MathNet.Numerics.Distributions
/// <param name="rng">The random number generator to use.</param>
/// <param name="mean">The mean of the normal distribution from which to generate samples.</param>
/// <param name="stddev">The standard deviation of the normal distribution from which to generate samples.</param>
/// <returns>a sequence of samples from the distribution.</returns>
public static IEnumerable<double> 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 <i>Box-Muller</i> algorithm.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="r2">The second random number.</param>
/// <param name="r2">A second random number from the standard normal distribution computed as a side product.</param>
/// <returns>a random number from the standard normal distribution.</returns>
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;

91
src/Managed/Distributions/Discrete/Bernoulli.cs
View File

@ -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<int> Samples() { throw new NotImplementedException(); }
#endregion
public static int Sample(System.Random rng, double p) { throw new NotImplementedException(); }
public static IEnumerable<int> Samples(System.Random rng, double p) { throw new NotImplementedException(); }
}
}
*/
// <copyright file="Bernoulli.cs" company="Math.NET">
// 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.
// </copyright>

38
src/Managed/Distributions/IContinuousDistribution.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="IContinuousDistribution.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
@ -35,14 +35,50 @@ namespace MathNet.Numerics.Distributions
/// </summary>
public interface IContinuousDistribution : IDistribution
{
/// <summary>
/// Gets the mode of the distribution.
/// </summary>
double Mode { get; }
/// <summary>
/// Gets the median of the distribution.
/// </summary>
double Median { get; }
/// <summary>
/// Gets the smallest element in the domain of the distributions which can be represented by a double.
/// </summary>
double Minimum { get; }
/// <summary>
/// Gets the largest element in the domain of the distributions which can be represented by a double.
/// </summary>
double Maximum { get; }
/// <summary>
/// The probability density of the distribution.
/// </summary>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
double Density(double x);
/// <summary>
/// The log probability density of the distribution.
/// </summary>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
double DensityLn(double x);
/// <summary>
/// Draws a random sample from the distribution.
/// </summary>
/// <returns>a sample from the distribution.</returns>
double Sample();
/// <summary>
/// Draws a sequence of random samples from the distribution.
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
IEnumerable<double> Samples();
}
}

14
src/Managed/Distributions/IDiscreteDistribution.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="IDiscreteDistribution.cs" company="Math.NET">
// 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
{
/// <summary>
/// The mode of the distribution.
/// Gets the mode of the distribution.
/// </summary>
int Mode { get; }
/// <summary>
/// The median of the distribution.
/// Gets the median of the distribution.
/// </summary>
int Median { get; }
/// <summary>
/// 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.
/// </summary>
int Minimum { get; }
/// <summary>
/// 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.
/// </summary>
int Maximum { get; }
@ -59,22 +59,26 @@ namespace MathNet.Numerics.Distributions
/// Computes values of the probability mass function.
/// </summary>
/// <param name="k">The location in the domain where we want to evaluate the probability mass function.</param>
/// <returns>the probability mass at location <paramref name="k"/>.</returns>
double Probability(int k);
/// <summary>
/// Computes values of the log probability mass function.
/// </summary>
/// <param name="k">The location in the domain where we want to evaluate the log probability mass function.</param>
/// <returns>the log probability mass at location <paramref name="k"/>.</returns>
double ProbabilityLn(int k);
/// <summary>
/// Draws a random sample from the distribution.
/// </summary>
/// <returns>a sample from the distribution.</returns>
int Sample();
/// <summary>
/// Draws a sequence of random samples from the distribution.
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
IEnumerable<int> Samples();
}
}

14
src/Managed/Distributions/IDistribution.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="IDistribution.cs" company="Math.NET">
// 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; }
/// <summary>
/// The mean of the distribution.
/// Gets the mean of the distribution.
/// </summary>
double Mean { get; }
/// <summary>
/// The variance of the distribution.
/// Gets the variance of the distribution.
/// </summary>
double Variance { get; }
/// <summary>
/// The standard deviation of the distribution.
/// Gets the standard deviation of the distribution.
/// </summary>
double StdDev { get; }
/// <summary>
/// The entropy of the distribution.
/// Gets the entropy of the distribution.
/// </summary>
double Entropy { get; }
/// <summary>
/// The skewness of the distribution.
/// Gets the skewness of the distribution.
/// </summary>
double Skewness { get; }
/// <summary>
/// Computes the cumulative distribution function (cdf) for this probability distribution.
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
double CumulativeDistribution(double x);
}
}

2
src/Managed/SiConstants.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="SiConstants.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//

2
src/Managed/SiPrefixes.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="SiPrefixes.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//

18
src/Managed/SpecialFunctions.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="SpecialFunctions.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
@ -34,21 +34,5 @@ namespace MathNet.Numerics
/// </summary>
public static partial class SpecialFunctions
{
/// <summary>
/// A helper function to evaluate polynomials fast.
/// </summary>
/// <param name="poly">The coefficients of the polynomial.</param>
/// <param name="z">The location where to evaluate the polynomial at.</param>
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;
}
}
}

222
src/Managed/SpecialFunctions/Erf.cs
View File

@ -1,4 +1,4 @@
// <copyright file="Combinatorics.cs" company="Math.NET">
// <copyright file="Erf.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
@ -30,11 +30,14 @@ namespace MathNet.Numerics
{
using System;
/// <summary>
/// This partial implementation of the SpecialFunctions class contains all methods related to the error function.
/// </summary>
public partial class SpecialFunctions
{
/// <summary>Calculates the error function.</summary>
/// <param name="x">The value to evaluate.</param>
/// <returns>The error function evaluated at given value.</returns>
/// <returns>the error function evaluated at given value.</returns>
/// <remarks>
/// <list type="bullet">
/// <item>returns 1 if <c>x == Double.PositiveInfinity</c>.</item>
@ -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
/// <summary>Calculates the complementary error function.</summary>
/// <param name="x">The value to evaluate.</param>
/// <returns>The complementary error function evaluated at given value.</returns>
/// <returns>the complementary error function evaluated at given value.</returns>
/// <remarks>
/// <list type="bullet">
/// <item>returns 0 if <c>x == Double.PositiveInfinity</c>.</item>
@ -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);
}
/// <summary>
/// Implementation of the error function.
/// </summary>
/// <param name="z">Where to evaluate the error function.</param>
/// <param name="invert">Whether to compute 1 - the error function.</param>
/// <returns>the error function.</returns>
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;
}
///<summary>Calculates the complementary inverse error function evaluated at z.</summary>
/// <returns>The complementary inverse error function evaluated at given value.</returns>
/// <remarks> We have tested this implementation against the arbitrary precision mpmath library
/// and found cases where we can only guarantee 9 significant figures correct.
/// <list type="bullet">
/// <item>returns Double.PositiveInfinity if <c>z &lt;= 0.0</c>.</item>
/// <item>returns Double.NegativeInfinity if <c>z &gt;= 2.0</c>.</item>
/// </list>
/// </remarks>
///<summary>Calculates the complementary inverse error function evaluated at z.</summary>
///<param name="z">value to evaluate.</param>
///<returns>the complementary inverse error function evaluated at Z.</returns>
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);
}
/// <summary>
/// The implementation of the inverse error function.
/// </summary>
/// <param name="p">First intermediate parameter.</param>
/// <param name="q">Second intermediate parameter.</param>
/// <param name="s">Third intermediate parameter.</param>
/// <returns>the inverse error function.</returns>
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;
}
/// <summary>
/// A helper function to evaluate polynomials fast.
/// </summary>
/// <param name="poly">The coefficients of the polynomial.</param>
/// <param name="z">The location where to evaluate the polynomial at.</param>
/// <returns>the evaluation of the polynomial.</returns>
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;
}
}
}

6
src/Native.UnitTests/Native.UnitTests.csproj
View File

@ -57,6 +57,9 @@
<Reference Include="System.Xml" />
</ItemGroup>
<ItemGroup>
<Compile Include="..\Managed.UnitTests\AssertHelpers.cs">
<Link>AssertHelpers.cs</Link>
</Compile>
<Compile Include="..\Managed.UnitTests\CombinatoricsTests\CombinatoricsCountingTest.cs">
<Link>CombinatoricsTests\CombinatoricsCountingTest.cs</Link>
</Compile>
@ -66,6 +69,9 @@
<Compile Include="..\Managed.UnitTests\DistributionTests\Continuous\NormalTests.cs">
<Link>DistributionTests\Continuous\NormalTests.cs</Link>
</Compile>
<Compile Include="..\Managed.UnitTests\PrecisionTest.cs">
<Link>PrecisionTest.cs</Link>
</Compile>
<Compile Include="..\Managed.UnitTests\SpecialFunctionsTest\ErfTests.cs">
<Link>SpecialFunctionsTest\ErfTests.cs</Link>
</Compile>

Write Preview
Loading…
Cancel
Save