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

Distributions: adapt Chi, ChiSquared

optimization-1
Christoph Ruegg 13 years ago
parent
221da13f42
commit
8075fd7f27
4 changed files with 147 additions and 73 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. 91
      src/Numerics/Distributions/Chi.cs
  2. 99
      src/Numerics/Distributions/ChiSquared.cs
  3. 12
      src/UnitTests/DistributionTests/Continuous/ChiSquareTests.cs
  4. 18
      src/UnitTests/DistributionTests/Continuous/ChiTests.cs

91
src/Numerics/Distributions/Chi.cs
View File

@ -82,16 +82,6 @@ namespace MathNet.Numerics.Distributions
return "Chi(k = " + _freedom + ")";
}
/// <summary>
/// Checks whether the parameters of the distribution are valid.
/// </summary>
/// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
static bool IsValidParameterSet(double freedom)
{
return freedom > 0.0;
}
/// <summary>
/// Sets the parameters of the distribution after checking their validity.
/// </summary>
@ -99,7 +89,7 @@ namespace MathNet.Numerics.Distributions
/// <exception cref="ArgumentOutOfRangeException">When the parameters are out of range.</exception>
void SetParameters(double freedom)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(freedom))
if (freedom <= 0.0 || Double.IsNaN(freedom))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
@ -214,6 +204,7 @@ namespace MathNet.Numerics.Distributions
/// </summary>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="PDF"/>
public double Density(double x)
{
return (Math.Pow(2.0, 1.0 - (_freedom/2.0))*Math.Pow(x, _freedom - 1.0)*Math.Exp(-x*x/2.0))/SpecialFunctions.Gamma(_freedom/2.0);
@ -224,6 +215,7 @@ namespace MathNet.Numerics.Distributions
/// </summary>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="PDFLn"/>
public double DensityLn(double x)
{
return ((1.0 - (_freedom/2.0))*Math.Log(2.0)) + ((_freedom - 1.0)*Math.Log(x)) - (x*x/2.0) - SpecialFunctions.GammaLn(_freedom/2.0);
@ -234,11 +226,34 @@ namespace MathNet.Numerics.Distributions
/// </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>
/// <seealso cref="CDF"/>
public double CumulativeDistribution(double x)
{
return SpecialFunctions.GammaLowerIncomplete(_freedom/2.0, x*x/2.0)/SpecialFunctions.Gamma(_freedom/2.0);
}
/// <summary>
/// Generates a sample from the Chi distribution.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public double Sample()
{
return SampleUnchecked(_random, (int) _freedom);
}
/// <summary>
/// Generates a sequence of samples from the Chi distribution.
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples()
{
var freedom = (int)_freedom;
while (true)
{
yield return SampleUnchecked(_random, freedom);
}
}
/// <summary>
/// Samples the distribution.
/// </summary>
@ -257,25 +272,45 @@ namespace MathNet.Numerics.Distributions
}
/// <summary>
/// Generates a sample from the Chi distribution.
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public double Sample()
/// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="Density"/>
public static double PDF(double freedom, double x)
{
return SampleUnchecked(_random, (int) _freedom);
if (freedom <= 0.0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
return (Math.Pow(2.0, 1.0 - (freedom/2.0))*Math.Pow(x, freedom - 1.0)*Math.Exp(-x*x/2.0))/SpecialFunctions.Gamma(freedom/2.0);
}
/// <summary>
/// Generates a sequence of samples from the Chi distribution.
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples()
/// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="DensityLn"/>
public static double PDFLn(double freedom, double x)
{
var freedom = (int)_freedom;
while (true)
{
yield return SampleUnchecked(_random, freedom);
}
if (freedom <= 0.0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
return ((1.0 - (freedom/2.0))*Math.Log(2.0)) + ((freedom - 1.0)*Math.Log(x)) - (x*x/2.0) - SpecialFunctions.GammaLn(freedom/2.0);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
public static double CDF(double freedom, double x)
{
if (freedom <= 0.0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
return SpecialFunctions.GammaLowerIncomplete(freedom/2.0, x*x/2.0)/SpecialFunctions.Gamma(freedom/2.0);
}
/// <summary>
@ -286,10 +321,7 @@ namespace MathNet.Numerics.Distributions
/// <returns>a sample from the distribution.</returns>
public static double Sample(System.Random rnd, int freedom)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(freedom))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
if (freedom <= 0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
return SampleUnchecked(rnd, freedom);
}
@ -302,10 +334,7 @@ namespace MathNet.Numerics.Distributions
/// <returns>a sequence of samples from the distribution.</returns>
public static IEnumerable<double> Samples(System.Random rnd, int freedom)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(freedom))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
if (freedom <= 0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
while (true)
{

99
src/Numerics/Distributions/ChiSquared.cs
View File

@ -80,16 +80,6 @@ namespace MathNet.Numerics.Distributions
return "ChiSquared(k = " + _freedom + ")";
}
/// <summary>
/// Checks whether the parameters of the distribution are valid.
/// </summary>
/// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
/// <returns><c>true</c> when the parameters are valid, <c>false</c> otherwise.</returns>
static bool IsValidParameterSet(double freedom)
{
return freedom > 0;
}
/// <summary>
/// Sets the parameters of the distribution after checking their validity.
/// </summary>
@ -97,7 +87,7 @@ namespace MathNet.Numerics.Distributions
/// <exception cref="ArgumentOutOfRangeException">When the parameters are out of range.</exception>
void SetParameters(double freedom)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(freedom))
if (freedom <= 0.0 || Double.IsNaN(freedom))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
@ -200,9 +190,10 @@ namespace MathNet.Numerics.Distributions
/// </summary>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="PDF"/>
public double Density(double x)
{
return (Math.Pow(x, (_freedom / 2.0) - 1.0) * Math.Exp(-x / 2.0)) / (Math.Pow(2.0, _freedom / 2.0) * SpecialFunctions.Gamma(_freedom / 2.0));
return (Math.Pow(x, (_freedom/2.0) - 1.0)*Math.Exp(-x/2.0))/(Math.Pow(2.0, _freedom/2.0)*SpecialFunctions.Gamma(_freedom/2.0));
}
/// <summary>
@ -210,9 +201,10 @@ namespace MathNet.Numerics.Distributions
/// </summary>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="PDFLn"/>
public double DensityLn(double x)
{
return (-x / 2.0) + (((_freedom / 2.0) - 1.0) * Math.Log(x)) - ((_freedom / 2.0) * Math.Log(2)) - SpecialFunctions.GammaLn(_freedom / 2.0);
return (-x/2.0) + (((_freedom/2.0) - 1.0)*Math.Log(x)) - ((_freedom/2.0)*Math.Log(2)) - SpecialFunctions.GammaLn(_freedom/2.0);
}
/// <summary>
@ -220,9 +212,31 @@ namespace MathNet.Numerics.Distributions
/// </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>
/// <seealso cref="CDF"/>
public double CumulativeDistribution(double x)
{
return SpecialFunctions.GammaLowerIncomplete(_freedom / 2.0, x / 2.0) / SpecialFunctions.Gamma(_freedom / 2.0);
return SpecialFunctions.GammaLowerIncomplete(_freedom/2.0, x/2.0)/SpecialFunctions.Gamma(_freedom/2.0);
}
/// <summary>
/// Generates a sample from the <c>ChiSquare</c> distribution.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public double Sample()
{
return SampleUnchecked(_random, _freedom);
}
/// <summary>
/// Generates a sequence of samples from the <c>ChiSquare</c> distribution.
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples()
{
while (true)
{
yield return SampleUnchecked(_random, _freedom);
}
}
/// <summary>
@ -237,7 +251,7 @@ namespace MathNet.Numerics.Distributions
if (Math.Floor(freedom) == freedom && freedom < Int32.MaxValue)
{
double sum = 0;
var n = (int) freedom;
var n = (int)freedom;
for (var i = 0; i < n; i++)
{
sum += Math.Pow(Normal.Sample(rnd, 0.0, 1.0), 2);
@ -247,28 +261,49 @@ namespace MathNet.Numerics.Distributions
//Call the gamma function (see http://en.wikipedia.org/wiki/Gamma_distribution#Specializations
//for a justification)
return Gamma.SampleUnchecked(rnd, freedom/2.0, .5);
return Gamma.SampleUnchecked(rnd, freedom / 2.0, .5);
}
/// <summary>
/// Generates a sample from the <c>ChiSquare</c> distribution.
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public double Sample()
/// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="Density"/>
public static double PDF(double freedom, double x)
{
return SampleUnchecked(_random, _freedom);
if (freedom <= 0.0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
return (Math.Pow(x, (freedom/2.0) - 1.0)*Math.Exp(-x/2.0))/(Math.Pow(2.0, freedom/2.0)*SpecialFunctions.Gamma(freedom/2.0));
}
/// <summary>
/// Generates a sequence of samples from the <c>ChiSquare</c> distribution.
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples()
/// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="DensityLn"/>
public static double PDFLn(double freedom, double x)
{
while (true)
{
yield return SampleUnchecked(_random, _freedom);
}
if (freedom <= 0.0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
return (-x/2.0) + (((freedom/2.0) - 1.0)*Math.Log(x)) - ((freedom/2.0)*Math.Log(2)) - SpecialFunctions.GammaLn(freedom/2.0);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
public static double CDF(double freedom, double x)
{
if (freedom <= 0.0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
return SpecialFunctions.GammaLowerIncomplete(freedom/2.0, x/2.0)/SpecialFunctions.Gamma(freedom/2.0);
}
/// <summary>
@ -279,10 +314,7 @@ namespace MathNet.Numerics.Distributions
/// <returns>a sample from the distribution. </returns>
public static double Sample(System.Random rnd, double freedom)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(freedom))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
if (freedom <= 0.0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
return SampleUnchecked(rnd, freedom);
}
@ -295,10 +327,7 @@ namespace MathNet.Numerics.Distributions
/// <returns>a sample from the distribution. </returns>
public static IEnumerable<double> Samples(System.Random rnd, double freedom)
{
if (Control.CheckDistributionParameters && !IsValidParameterSet(freedom))
{
throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
}
if (freedom <= 0.0) throw new ArgumentOutOfRangeException("freedom", Resources.InvalidDistributionParameters);
while (true)
{

12
src/UnitTests/DistributionTests/Continuous/ChiSquareTests.cs
View File

@ -241,7 +241,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
public void ValidateDensity(double dof, double x)
{
var n = new ChiSquared(dof);
Assert.AreEqual((Math.Pow(x, (dof / 2.0) - 1.0) * Math.Exp(-x / 2.0)) / (Math.Pow(2.0, dof / 2.0) * SpecialFunctions.Gamma(dof / 2.0)), n.Density(x));
double expected = (Math.Pow(x, (dof / 2.0) - 1.0) * Math.Exp(-x / 2.0)) / (Math.Pow(2.0, dof / 2.0) * SpecialFunctions.Gamma(dof / 2.0));
Assert.AreEqual(expected, n.Density(x));
Assert.AreEqual(expected, ChiSquared.PDF(dof, x));
}
/// <summary>
@ -276,7 +278,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
public void ValidateDensityLn(double dof, double x)
{
var n = new ChiSquared(dof);
Assert.AreEqual((-x / 2.0) + (((dof / 2.0) - 1.0) * Math.Log(x)) - ((dof / 2.0) * Math.Log(2)) - SpecialFunctions.GammaLn(dof / 2.0), n.DensityLn(x));
double expected = (-x / 2.0) + (((dof / 2.0) - 1.0) * Math.Log(x)) - ((dof / 2.0) * Math.Log(2)) - SpecialFunctions.GammaLn(dof / 2.0);
Assert.AreEqual(expected, n.DensityLn(x));
Assert.AreEqual(expected, ChiSquared.PDFLn(dof, x));
}
/// <summary>
@ -350,7 +354,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
public void ValidateCumulativeDistribution(double dof, double x)
{
var n = new ChiSquared(dof);
Assert.AreEqual(SpecialFunctions.GammaLowerIncomplete(dof / 2.0, x / 2.0) / SpecialFunctions.Gamma(dof / 2.0), n.CumulativeDistribution(x));
double expected = SpecialFunctions.GammaLowerIncomplete(dof / 2.0, x / 2.0) / SpecialFunctions.Gamma(dof / 2.0);
Assert.AreEqual(expected, n.CumulativeDistribution(x));
Assert.AreEqual(expected, ChiSquared.CDF(dof, x));
}
}
}

18
src/UnitTests/DistributionTests/Continuous/ChiTests.cs
View File

@ -3,7 +3,9 @@
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
// Copyright (c) 2009-2010 Math.NET
//
// Copyright (c) 2009-2013 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
@ -12,8 +14,10 @@
// 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
@ -237,7 +241,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
public void ValidateDensity(double dof, double x)
{
var n = new Chi(dof);
Assert.AreEqual((Math.Pow(2.0, 1.0 - (dof / 2.0)) * Math.Pow(x, dof - 1.0) * Math.Exp(-x * (x / 2.0))) / SpecialFunctions.Gamma(dof / 2.0), n.Density(x));
double expected = (Math.Pow(2.0, 1.0 - (dof / 2.0)) * Math.Pow(x, dof - 1.0) * Math.Exp(-x * (x / 2.0))) / SpecialFunctions.Gamma(dof / 2.0);
Assert.AreEqual(expected, n.Density(x));
Assert.AreEqual(expected, Chi.PDF(dof, x));
}
/// <summary>
@ -272,7 +278,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
public void ValidateDensityLn(double dof, double x)
{
var n = new Chi(dof);
Assert.AreEqual(((1.0 - (dof / 2.0)) * Math.Log(2.0)) + ((dof - 1.0) * Math.Log(x)) - (x * (x / 2.0)) - SpecialFunctions.GammaLn(dof / 2.0), n.DensityLn(x));
double expected = ((1.0 - (dof / 2.0)) * Math.Log(2.0)) + ((dof - 1.0) * Math.Log(x)) - (x * (x / 2.0)) - SpecialFunctions.GammaLn(dof / 2.0);
Assert.AreEqual(expected, n.DensityLn(x));
Assert.AreEqual(expected, Chi.PDFLn(dof, x));
}
/// <summary>
@ -328,7 +336,9 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
public void ValidateCumulativeDistribution(double dof, double x)
{
var n = new Chi(dof);
Assert.AreEqual(SpecialFunctions.GammaLowerIncomplete(dof / 2.0, x * x / 2.0) / SpecialFunctions.Gamma(dof / 2.0), n.CumulativeDistribution(x));
double expected = SpecialFunctions.GammaLowerIncomplete(dof / 2.0, x * x / 2.0) / SpecialFunctions.Gamma(dof / 2.0);
Assert.AreEqual(expected, n.CumulativeDistribution(x));
Assert.AreEqual(expected, Chi.CDF(dof, x));
}
}
}

Write Preview
Loading…
Cancel
Save