From dee7fc6cb9de32cdd0b23f2cc203f3ea7846afaa Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 1 May 2014 20:34:14 +0200 Subject: [PATCH] =?UTF-8?q?add=20Fr=C3=A9chet=20distribution?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- doc/source/univariate.rst | 20 ++++++++++ src/Distributions.jl | 1 + src/univariate/frechet.jl | 81 +++++++++++++++++++++++++++++++++++++++ src/univariates.jl | 1 + test/univariate.jl | 12 ++++++ 5 files changed, 115 insertions(+) create mode 100644 src/univariate/frechet.jl diff --git a/doc/source/univariate.rst b/doc/source/univariate.rst index 54def0dab..af4a059d1 100644 --- a/doc/source/univariate.rst +++ b/doc/source/univariate.rst @@ -205,6 +205,7 @@ List of Distributions - :ref:`erlang` - :ref:`exponential` - :ref:`fdist` + - :ref:`frechet` - :ref:`gamma` - :ref:`gumbel` - :ref:`inversegamma` @@ -494,6 +495,25 @@ The probability density function of an `F distribution 0 and scale θ>0 is + +.. math:: + + f(x; k, \theta) = \frac{k}{\theta} \left( \frac{x}{\theta} \right)^{-k-1} e^{-(x/\theta)^{-k}}, + \quad x > 0 + +.. code-block:: julia + + Frechet(k) # Fréchet distribution with shape k and unit scale + Frechet(k, s) # Fréchet distribution with shape k and scale s + + .. _gamma: Gamma Distribution diff --git a/src/Distributions.jl b/src/Distributions.jl index ef0bd4fa6..97cfa2cad 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -60,6 +60,7 @@ export Erlang, Exponential, FDist, + Frechet, Gamma, GenericMvNormal, GenericMvNormalCanon, diff --git a/src/univariate/frechet.jl b/src/univariate/frechet.jl new file mode 100644 index 000000000..ed277f44c --- /dev/null +++ b/src/univariate/frechet.jl @@ -0,0 +1,81 @@ +# Fréchet Distribution + +immutable Frechet <: ContinuousUnivariateDistribution + shape::Float64 + scale::Float64 + function Frechet(sh::Real, sc::Real) + zero(sh) < sh && zero(sc) < sc || error("Both shape and scale must be positive") + new(float64(sh), float64(sc)) + end +end + +Frechet(sh::Real) = Frechet(sh, 1.0) + +## Support +@continuous_distr_support Frechet 0.0 Inf + +## Properties +mean(d::Frechet) = d.shape > 1.0 ? d.scale * gamma(1.0 - 1.0 / d.shape) : Inf +median(d::Frechet) = d.scale * log(2)^(-1.0 / d.shape) + +mode(d::Frechet) = (ik = -1.0/d.shape; d.scale * (1.0-ik)^ik) + +var(d::Frechet) = d.shape > 2.0 ? d.scale^2 * gamma(1.0 - 2.0 / d.shape) - mean(d)^2 : NaN + +function skewness(d::Frechet) + d.shape <= 3.0 && return NaN + tmp_mean = mean(d) + tmp_var = var(d) + tmp = gamma(1.0 - 3.0 / d.shape) * d.scale^3 + tmp -= 3.0 * tmp_mean * tmp_var + tmp -= tmp_mean^3 + return tmp / tmp_var / sqrt(tmp_var) +end + +function kurtosis(d::Frechet) + d.shape <= 4.0 && return NaN + λ, k = d.scale, d.shape + μ = mean(d) + σ = std(d) + γ = skewness(d) + den = λ^4 * gamma(1.0 - 4.0 / k) - + 4.0 * γ * σ^3 * μ - + 6.0 * μ^2 * σ^2 - μ^4 + num = σ^4 + return den / num - 3.0 +end + +function entropy(d::Frechet) + λ, k = d.scale, d.shape + return (k + 1.0) * (log(λ) - digamma(1.0)/k) - log(λ * k) + 1.0 +end + +## Functions +function pdf(d::Frechet, x::Real) + x < zero(x) && return 0.0 + a = d.scale/x + d.shape/d.scale * a^(d.shape+1.0) * exp(-a^d.shape) +end +function logpdf(d::Frechet, x::Real) + x < zero(x) && return -Inf + a = d.scale/x + log(d.shape/d.scale) + (d.shape+1.0)*log(a) - a^d.shape +end + +cdf(d::Frechet, x::Real) = x <= zero(x) ? 0.0 : exp(-((d.scale / x)^d.shape)) +ccdf(d::Frechet, x::Real) = x <= zero(x) ? 1.0 : -expm1(-((d.scale / x)^d.shape)) +logcdf(d::Frechet, x::Real) = x <= zero(x) ? -Inf : -(d.scale / x)^d.shape +logccdf(d::Frechet, x::Real) = x <= zero(x) ? 0.0 : log1mexp(-((d.scale / x)^d.shape)) + +quantile(d::Frechet, p::Real) = @checkquantile p d.scale*(-log(p))^(-1/d.shape) +cquantile(d::Frechet, p::Real) = @checkquantile p d.scale*(-log1p(-p))^(-1/d.shape) +invlogcdf(d::Frechet, lp::Real) = lp > zero(lp) ? NaN : d.scale*(-lp)^(-1/d.shape) +invlogccdf(d::Frechet, lp::Real) = lp > zero(lp) ? NaN : d.scale*(-log1mexp(lp))^(-1/d.shape) + + +function gradlogpdf(d::Frechet, x::Real) + insupport(Frechet, x) ? -(d.shape + 1.0) / x + d.shape * (d.scale^d.shape) * x^(-d.shape - 1.0) : 0.0 +end + +## Sampling +rand(d::Frechet) = d.scale*Base.Random.randmtzig_exprnd()^(-1/d.shape) diff --git a/src/univariates.jl b/src/univariates.jl index eafc8d66a..4c282fa00 100644 --- a/src/univariates.jl +++ b/src/univariates.jl @@ -123,6 +123,7 @@ for finame in ["arcsine.jl", "empirical.jl", "exponential.jl", "fdist.jl", + "frechet.jl", "gamma.jl", "edgeworth.jl", "erlang.jl", diff --git a/test/univariate.jl b/test/univariate.jl index 2298d23e0..f278a59f8 100644 --- a/test/univariate.jl +++ b/test/univariate.jl @@ -44,6 +44,18 @@ distlist = [Arcsine(), FDist(9, 9), FDist(9, 21), FDist(21, 9), + Frechet(0.23,0.1), + Frechet(2.3,0.1), + Frechet(23.0,0.1), + Frechet(230.0,0.1), + Frechet(0.23), + Frechet(2.3), + Frechet(23.0), + Frechet(230.0), + Frechet(0.23,10.0), + Frechet(2.3,10.0), + Frechet(23.0,10.0), + Frechet(230.0,10.0), Gamma(3.0, 2.0), Gamma(2.0, 3.0), Gamma(3.0, 3.0),