From b4592acd0297dde2d4cea340b1c368ca2e3965e0 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Mon, 30 Dec 2013 21:52:00 +0100 Subject: [PATCH] Distributions: static PDF/CDF functions for StudentT, Stable, Weibull #188 --- src/Numerics/Distributions/Stable.cs | 188 +++++++++--------- src/Numerics/Distributions/StudentT.cs | 145 ++++++++------ src/Numerics/Distributions/Weibull.cs | 107 +++++++--- .../Continuous/StudentTTests.cs | 19 -- 4 files changed, 256 insertions(+), 203 deletions(-) diff --git a/src/Numerics/Distributions/Stable.cs b/src/Numerics/Distributions/Stable.cs index c4897cd1..bf3b7013 100644 --- a/src/Numerics/Distributions/Stable.cs +++ b/src/Numerics/Distributions/Stable.cs @@ -37,10 +37,10 @@ namespace MathNet.Numerics.Distributions { /// /// Continuous Univariate Stable distribution. - /// A random variable is said to be stable (or to have a stable distribution) if it has - /// the property that a linear combination of two independent copies of the variable has + /// A random variable is said to be stable (or to have a stable distribution) if it has + /// the property that a linear combination of two independent copies of the variable has /// the same distribution, up to location and scale parameters. - /// For details about this distribution, see + /// For details about this distribution, see /// Wikipedia - Stable distribution. /// /// The distribution will use the by default.` @@ -58,7 +58,7 @@ namespace MathNet.Numerics.Distributions double _location; /// - /// Initializes a new instance of the class. + /// Initializes a new instance of the class. /// /// The stability (α) of the distribution. Range: 2 ≥ α > 0. /// The skewness (β) of the distribution. Range: 1 ≥ β ≥ -1. @@ -71,7 +71,7 @@ namespace MathNet.Numerics.Distributions } /// - /// Initializes a new instance of the class. + /// Initializes a new instance of the class. /// /// The stability (α) of the distribution. Range: 2 ≥ α > 0. /// The skewness (β) of the distribution. Range: 1 ≥ β ≥ -1. @@ -93,19 +93,6 @@ namespace MathNet.Numerics.Distributions return "Stable(α = " + _alpha + ", β = " + _beta + ", c = " + _scale + ", μ = " + _location + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The stability (α) of the distribution. Range: 2 ≥ α > 0. - /// The skewness (β) of the distribution. Range: 1 ≥ β ≥ -1. - /// The scale (c) of the distribution. Range: c > 0. - /// The location (μ) of the distribution. - /// true when the parameters are valid, false otherwise. - static bool IsValidParameterSet(double alpha, double beta, double scale, double location) - { - return alpha > 0.0 && alpha <= 2.0 && beta >= -1.0 && beta <= 1.0 && scale > 0.0 && !Double.IsNaN(location); - } - /// /// Sets the parameters of the distribution after checking their validity. /// @@ -116,7 +103,8 @@ namespace MathNet.Numerics.Distributions /// When the parameters are out of range. void SetParameters(double alpha, double beta, double scale, double location) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(alpha, beta, scale, location)) + if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0 + || Double.IsNaN(alpha) || Double.IsNaN(beta) || Double.IsNaN(scale) || Double.IsNaN(location)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } @@ -179,7 +167,7 @@ namespace MathNet.Numerics.Distributions { get { - if (_alpha <= 1) + if (_alpha <= 1d) { throw new NotSupportedException(); } @@ -195,7 +183,7 @@ namespace MathNet.Numerics.Distributions { get { - if (_alpha == 2) + if (_alpha == 2d) { return 2.0*_scale*_scale; } @@ -211,7 +199,7 @@ namespace MathNet.Numerics.Distributions { get { - if (_alpha == 2) + if (_alpha == 2d) { return Constants.Sqrt2*_scale; } @@ -237,7 +225,7 @@ namespace MathNet.Numerics.Distributions { get { - if (_alpha != 2) + if (_alpha != 2d) { throw new NotSupportedException(); } @@ -254,7 +242,7 @@ namespace MathNet.Numerics.Distributions { get { - if (_beta != 0) + if (_beta != 0d) { throw new NotSupportedException(); } @@ -311,40 +299,7 @@ namespace MathNet.Numerics.Distributions /// the density at . public double Density(double x) { - if (_alpha == 2) - { - return (new Normal(_location, StdDev)).Density(x); - } - - if (_alpha == 1 && _beta == 0) - { - return (new Cauchy(_location, _scale)).Density(x); - } - - if (_alpha == 0.5 && _beta == 1) - { - return LevyDensity(_scale, _location, x); - } - - throw new NotSupportedException(); - } - - /// - /// Computes the density of the Levy distribution. - /// - /// The scale (c) of the distribution. - /// The location (μ) of the distribution. - /// The location at which to compute the density. - /// the density at . - static double LevyDensity(double scale, double location, double x) - { - // The parameters scale and location must be correct - if (x < location) - { - throw new NotSupportedException(); - } - - return (Math.Sqrt(scale/Constants.Pi2)*Math.Exp(-scale/(2*(x - location))))/Math.Pow(x - location, 1.5); + return PDF(_alpha, _beta, _scale, _location, x); } /// @@ -354,7 +309,7 @@ namespace MathNet.Numerics.Distributions /// the log density at . public double DensityLn(double x) { - return Math.Log(Density(x)); + return PDFLn(_alpha, _beta, _scale, _location, x); } /// @@ -365,35 +320,7 @@ namespace MathNet.Numerics.Distributions /// Throws a not supported exception if Alpha != 2, (Alpha != 1 and Beta !=0), or (Alpha != 0.5 and Beta != 1) public double CumulativeDistribution(double x) { - if (_alpha == 2) - { - return (new Normal(_location, StdDev)).CumulativeDistribution(x); - } - - if (_alpha == 1 && _beta == 0) - { - return (new Cauchy(_location, _scale)).CumulativeDistribution(x); - } - - if (_alpha == 0.5 && _beta == 1) - { - return LevyCumulativeDistribution(_scale, _location, x); - } - - throw new NotSupportedException(); - } - - /// - /// Computes the cumulative distribution function of the Levy distribution. - /// - /// The scale (c) of the distribution. - /// The location (μ) of the distribution. - /// The location at which to compute the cumulative density. - /// the cumulative density at . - static double LevyCumulativeDistribution(double scale, double location, double x) - { - // The parameters scale and location must be correct - return SpecialFunctions.Erfc(Math.Sqrt(scale/(2*(x - location)))); + return CDF(_alpha, _beta, _scale, _location, x); } /// @@ -453,6 +380,83 @@ namespace MathNet.Numerics.Distributions } } + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The stability (α) of the distribution. Range: 2 ≥ α > 0. + /// The skewness (β) of the distribution. Range: 1 ≥ β ≥ -1. + /// The scale (c) of the distribution. Range: c > 0. + /// The location (μ) of the distribution. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double alpha, double beta, double scale, double location, double x) + { + if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0) + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + if (alpha == 2d) return Normal.PDF(location, Constants.Sqrt2*scale, x); + if (alpha == 1d && beta == 0d) return Cauchy.PDF(location, scale, x); + + if (alpha == 0.5d && beta == 1d && x >= location) + { + return (Math.Sqrt(scale/Constants.Pi2)*Math.Exp(-scale/(2*(x - location))))/Math.Pow(x - location, 1.5); + } + + throw new NotSupportedException(); + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The stability (α) of the distribution. Range: 2 ≥ α > 0. + /// The skewness (β) of the distribution. Range: 1 ≥ β ≥ -1. + /// The scale (c) of the distribution. Range: c > 0. + /// The location (μ) of the distribution. + /// The location at which to compute the density. + /// the log density at . + /// + public static double PDFLn(double alpha, double beta, double scale, double location, double x) + { + if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0) + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + if (alpha == 2d) return Normal.PDFLn(location, Constants.Sqrt2*scale, x); + if (alpha == 1d && beta == 0d) return Cauchy.PDFLn(location, scale, x); + + if (alpha == 0.5d && beta == 1d && x >= location) + { + return (Math.Log(scale/Constants.Pi2))/2 - scale/(2*(x - location)) - 1.5*Math.Log(x - location); + } + + throw new NotSupportedException(); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The location at which to compute the cumulative distribution function. + /// The stability (α) of the distribution. Range: 2 ≥ α > 0. + /// The skewness (β) of the distribution. Range: 1 ≥ β ≥ -1. + /// The scale (c) of the distribution. Range: c > 0. + /// The location (μ) of the distribution. + /// the cumulative distribution at location . + /// + public static double CDF(double alpha, double beta, double scale, double location, double x) + { + if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0) + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + if (alpha == 2d) return Normal.CDF(location, Constants.Sqrt2*scale, x); + if (alpha == 1d && beta == 0d) return Cauchy.CDF(location, scale, x); + + if (alpha == 0.5d && beta == 1d) + { + return SpecialFunctions.Erfc(Math.Sqrt(scale/(2*(x - location)))); + } + + throw new NotSupportedException(); + } /// /// Generates a sample from the distribution. @@ -465,10 +469,8 @@ namespace MathNet.Numerics.Distributions /// a sample from the distribution. public static double Sample(System.Random rnd, double alpha, double beta, double scale, double location) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(alpha, beta, scale, location)) - { + if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } return SampleUnchecked(rnd, alpha, beta, scale, location); } @@ -484,10 +486,8 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rnd, double alpha, double beta, double scale, double location) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale, scale, location)) - { + if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } while (true) { diff --git a/src/Numerics/Distributions/StudentT.cs b/src/Numerics/Distributions/StudentT.cs index 92bb1210..df0ec871 100644 --- a/src/Numerics/Distributions/StudentT.cs +++ b/src/Numerics/Distributions/StudentT.cs @@ -38,7 +38,7 @@ namespace MathNet.Numerics.Distributions /// /// Continuous Univariate Student's T-distribution. /// Implements the univariate Student t-distribution. For details about this - /// distribution, see + /// distribution, see /// /// Wikipedia - Student's t-distribution. /// @@ -50,7 +50,7 @@ namespace MathNet.Numerics.Distributions /// Gamma((dof+1)/2) (1 + (x - mu)^2 / (scale * scale * dof))^(-(dof+1)/2) / /// (Gamma(dof/2)*Sqrt(dof*pi*scale)). /// The distribution will use the by - /// default. Users can get/set the random number generator by using the + /// default. Users can get/set the random number generator by using the /// property. /// The statistics classes will check all the incoming parameters /// whether they are in the allowed range. This might involve heavy @@ -66,18 +66,17 @@ namespace MathNet.Numerics.Distributions /// /// Initializes a new instance of the StudentT class. This is a Student t-distribution with location 0.0 - /// scale 1.0 and degrees of freedom 1. The distribution will - /// be initialized with the default random number generator. + /// scale 1.0 and degrees of freedom 1. /// public StudentT() - : this(0.0, 1.0, 1.0) { + _random = MersenneTwister.Default; + SetParameters(0.0, 1.0, 1.0); } /// /// Initializes a new instance of the StudentT class with a particular location, scale and degrees of - /// freedom. The distribution will - /// be initialized with the default random number generator. + /// freedom. /// /// The location (μ) of the distribution. /// The scale (σ) of the distribution. Range: σ > 0. @@ -90,8 +89,7 @@ namespace MathNet.Numerics.Distributions /// /// Initializes a new instance of the StudentT class with a particular location, scale and degrees of - /// freedom. The distribution will - /// be initialized with the default random number generator. + /// freedom. /// /// The location (μ) of the distribution. /// The scale (σ) of the distribution. Range: σ > 0. @@ -112,18 +110,6 @@ namespace MathNet.Numerics.Distributions return "StudentT(μ = " + _location + ", σ = " + _scale + ", ν = " + _freedom + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The location (μ) of the distribution. - /// The scale (σ) of the distribution. Range: σ > 0. - /// The degrees of freedom (ν) for the distribution. Range: ν > 0. - /// true when the parameters are valid, false otherwise. - static bool IsValidParameterSet(double location, double scale, double freedom) - { - return scale > 0.0 && freedom > 0.0 && !Double.IsNaN(location); - } - /// /// Sets the parameters of the distribution after checking their validity. /// @@ -133,7 +119,7 @@ namespace MathNet.Numerics.Distributions /// When the parameters are out of range. void SetParameters(double location, double scale, double freedom) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale, freedom)) + if (scale <= 0.0 || freedom <= 0.0 || Double.IsNaN(scale) || Double.IsNaN(location) || Double.IsNaN(freedom)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } @@ -295,17 +281,7 @@ namespace MathNet.Numerics.Distributions /// the density at . public double Density(double x) { - // TODO JVG we can probably do a better job for Cauchy special case - if (_freedom >= 1e+8d) - { - return Normal.PDF(_location, _scale, x); - } - - var d = (x - _location)/_scale; - return Math.Exp(SpecialFunctions.GammaLn((_freedom + 1.0)/2.0) - SpecialFunctions.GammaLn(_freedom/2.0)) - *Math.Pow(1.0 + (d*d/_freedom), -0.5*(_freedom + 1.0)) - /Math.Sqrt(_freedom*Math.PI) - /_scale; + return PDF(_location, _scale, _freedom, x); } /// @@ -315,17 +291,7 @@ namespace MathNet.Numerics.Distributions /// the log density at . public double DensityLn(double x) { - // TODO JVG we can probably do a better job for Cauchy special case - if (_freedom >= 1e+8d) - { - return Normal.PDFLn(_location, _scale, x); - } - - var d = (x - _location)/_scale; - return SpecialFunctions.GammaLn((_freedom + 1.0)/2.0) - - (0.5*((_freedom + 1.0)*Math.Log(1.0 + (d*d/_freedom)))) - - SpecialFunctions.GammaLn(_freedom/2.0) - - (0.5*Math.Log(_freedom*Math.PI)) - Math.Log(_scale); + return PDFLn(_location, _scale, _freedom, x); } /// @@ -335,22 +301,13 @@ namespace MathNet.Numerics.Distributions /// the cumulative distribution at location . public double CumulativeDistribution(double x) { - // TODO JVG we can probably do a better job for Cauchy special case - if (Double.IsPositiveInfinity(_freedom)) - { - return Normal.CDF(_location, _scale, x); - } - - var k = (x - _location)/_scale; - var h = _freedom/(_freedom + (k*k)); - var ib = 0.5*SpecialFunctions.BetaRegularized(_freedom/2.0, 0.5, h); - return x <= _location ? ib : 1.0 - ib; + return CDF(_location, _scale, _freedom, x); } /// /// Samples student-t distributed random variables. /// - /// The algorithm is method 2 in section 5, chapter 9 + /// The algorithm is method 2 in section 5, chapter 9 /// in L. Devroye's "Non-Uniform Random Variate Generation" /// The random number generator to use. /// The location (μ) of the distribution. @@ -384,6 +341,74 @@ namespace MathNet.Numerics.Distributions } } + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The degrees of freedom (ν) for the distribution. Range: ν > 0. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double location, double scale, double freedom, double x) + { + if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + // TODO JVG we can probably do a better job for Cauchy special case + if (freedom >= 1e+8d) return Normal.PDF(location, scale, x); + + var d = (x - location)/scale; + return Math.Exp(SpecialFunctions.GammaLn((freedom + 1.0)/2.0) - SpecialFunctions.GammaLn(freedom/2.0)) + *Math.Pow(1.0 + (d*d/freedom), -0.5*(freedom + 1.0)) + /Math.Sqrt(freedom*Math.PI) + /scale; + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The degrees of freedom (ν) for the distribution. Range: ν > 0. + /// The location at which to compute the density. + /// the log density at . + /// + public static double PDFLn(double location, double scale, double freedom, double x) + { + if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + // TODO JVG we can probably do a better job for Cauchy special case + if (freedom >= 1e+8d) return Normal.PDFLn(location, scale, x); + + var d = (x - location)/scale; + return SpecialFunctions.GammaLn((freedom + 1.0)/2.0) + - (0.5*((freedom + 1.0)*Math.Log(1.0 + (d*d/freedom)))) + - SpecialFunctions.GammaLn(freedom/2.0) + - (0.5*Math.Log(freedom*Math.PI)) - Math.Log(scale); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The location at which to compute the cumulative distribution function. + /// The location (μ) of the distribution. + /// The scale (σ) of the distribution. Range: σ > 0. + /// The degrees of freedom (ν) for the distribution. Range: ν > 0. + /// the cumulative distribution at location . + /// + public static double CDF(double location, double scale, double freedom, double x) + { + if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + // TODO JVG we can probably do a better job for Cauchy special case + if (Double.IsPositiveInfinity(freedom)) return Normal.CDF(location, scale, x); + + var k = (x - location)/scale; + var h = freedom/(freedom + (k*k)); + var ib = 0.5*SpecialFunctions.BetaRegularized(freedom/2.0, 0.5, h); + return x <= location ? ib : 1.0 - ib; + } + /// /// Generates a sample from the Student t-distribution. /// @@ -394,10 +419,7 @@ namespace MathNet.Numerics.Distributions /// a sample from the distribution. public static double Sample(System.Random rnd, double location, double scale, double freedom) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale, freedom)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); return SampleUnchecked(rnd, location, scale, freedom); } @@ -412,10 +434,7 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rnd, double location, double scale, double freedom) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(location, scale, freedom)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); while (true) { diff --git a/src/Numerics/Distributions/Weibull.cs b/src/Numerics/Distributions/Weibull.cs index 8d282d85..81e559e2 100644 --- a/src/Numerics/Distributions/Weibull.cs +++ b/src/Numerics/Distributions/Weibull.cs @@ -4,7 +4,7 @@ // 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 @@ -37,12 +37,12 @@ namespace MathNet.Numerics.Distributions { /// /// Continuous Univariate Weibull distribution. - /// For details about this distribution, see + /// For details about this distribution, see /// Wikipedia - Weibull distribution. /// /// /// The Weibull distribution is parametrized by a shape and scale parameter. - /// The distribution will use the by default. + /// The distribution will use the by default. /// Users can get/set the random number generator by using the property. /// The statistics classes will check all the incoming parameters whether they are in the allowed /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters @@ -95,17 +95,6 @@ namespace MathNet.Numerics.Distributions return "Weibull(k = " + _shape + ", λ = " + _scale + ")"; } - /// - /// Checks whether the parameters of the distribution are valid. - /// - /// The shape (k) of the Weibull distribution. Range: k > 0. - /// The scale (λ) of the Weibull distribution. Range: λ > 0. - /// true when the parameters positive valid floating point numbers, false otherwise. - static bool IsValidParameterSet(double shape, double scale) - { - return shape > 0.0 && scale > 0.0; - } - /// /// Sets the parameters of the distribution after checking their validity. /// @@ -114,7 +103,7 @@ namespace MathNet.Numerics.Distributions /// When the parameters are out of range. void SetParameters(double shape, double scale) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) + if (shape <= 0.0 || scale <= 0.0 || Double.IsNaN(shape) || Double.IsNaN(scale)) { throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); } @@ -285,10 +274,7 @@ namespace MathNet.Numerics.Distributions /// the cumulative distribution at location . public double CumulativeDistribution(double x) { - if (x < 0.0) - { - return 0.0; - } + if (x < 0.0) return 0.0; return -SpecialFunctions.ExponentialMinusOne(-Math.Pow(x, _shape)*_scalePowShapeInv); } @@ -328,6 +314,79 @@ namespace MathNet.Numerics.Distributions } } + /// + /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x. + /// + /// The shape (k) of the Weibull distribution. Range: k > 0. + /// The scale (λ) of the Weibull distribution. Range: λ > 0. + /// The location at which to compute the density. + /// the density at . + /// + public static double PDF(double shape, double scale, double x) + { + if (shape <= 0.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + if (x >= 0.0) + { + if (x == 0.0 && shape == 1.0) + { + return shape/scale; + } + + return shape + *Math.Pow(x/scale, shape - 1.0) + *Math.Exp(-Math.Pow(x, shape)*Math.Pow(scale, -shape)) + /scale; + } + + return 0.0; + } + + /// + /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x). + /// + /// The shape (k) of the Weibull distribution. Range: k > 0. + /// The scale (λ) of the Weibull distribution. Range: λ > 0. + /// The location at which to compute the density. + /// the log density at . + /// + public static double PDFLn(double shape, double scale, double x) + { + if (shape <= 0.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + if (x >= 0.0) + { + if (x == 0.0 && shape == 1.0) + { + return Math.Log(shape) - Math.Log(scale); + } + + return Math.Log(shape) + + ((shape - 1.0)*Math.Log(x/scale)) + - (Math.Pow(x, shape)*Math.Pow(scale, -shape)) + - Math.Log(scale); + } + + return double.NegativeInfinity; + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The location at which to compute the cumulative distribution function. + /// The shape (k) of the Weibull distribution. Range: k > 0. + /// The scale (λ) of the Weibull distribution. Range: λ > 0. + /// the cumulative distribution at location . + /// + public static double CDF(double shape, double scale, double x) + { + if (shape <= 0.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + + if (x < 0.0) return 0.0; + + return -SpecialFunctions.ExponentialMinusOne(-Math.Pow(x, shape)*Math.Pow(scale, -shape)); + } + /// /// Generates a sample from the Weibull distribution. /// @@ -337,10 +396,7 @@ namespace MathNet.Numerics.Distributions /// a sample from the distribution. public static double Sample(System.Random rnd, double shape, double scale) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (shape <= 0.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); return SampleUnchecked(rnd, shape, scale); } @@ -354,10 +410,7 @@ namespace MathNet.Numerics.Distributions /// a sequence of samples from the distribution. public static IEnumerable Samples(System.Random rnd, double shape, double scale) { - if (Control.CheckDistributionParameters && !IsValidParameterSet(shape, scale)) - { - throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); - } + if (shape <= 0.0 || scale <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); while (true) { diff --git a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs index 2096fcdf..0f716cb6 100644 --- a/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs +++ b/src/UnitTests/DistributionTests/Continuous/StudentTTests.cs @@ -390,25 +390,6 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous ied.Take(5).ToArray(); } - /// - /// Fail sample static with bad parameters. - /// - [Test] - public void FailSampleStatic() - { - Assert.Throws(() => StudentT.Sample(new Random(0), Double.NaN, 1.0, Double.NaN)); - } - - /// - /// Fail sample sequence static with bad parameters. - /// - [Test] - public void FailSampleSequenceStatic() - { - var ied = StudentT.Samples(new Random(0), 0.0, 1.0, Double.NaN); - Assert.Throws(() => ied.Take(5).ToArray()); - } - /// /// Can sample. ///