Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updating sympy.stats.crv_types for extreme value distributions #16528

Merged
merged 17 commits into from Jun 16, 2019
Merged
2 changes: 1 addition & 1 deletion sympy/core/tests/test_args.py
Expand Up @@ -1359,7 +1359,7 @@ def test_sympy__stats__crv_types__GammaDistribution():

def test_sympy__stats__crv_types__GumbelDistribution():
from sympy.stats.crv_types import GumbelDistribution
assert _test_args(GumbelDistribution(1, 1))
assert _test_args(GumbelDistribution(1, 1, False))

def test_sympy__stats__crv_types__GompertzDistribution():
from sympy.stats.crv_types import GompertzDistribution
Expand Down
52 changes: 37 additions & 15 deletions sympy/stats/crv_types.py
Expand Up @@ -50,7 +50,7 @@
from sympy import (log, sqrt, pi, S, Dummy, Interval, sympify, gamma,
Piecewise, And, Eq, binomial, factorial, Sum, floor, Abs,
Lambda, Basic, lowergamma, erf, erfi, erfinv, I, hyper,
uppergamma, sinh, atan, Ne, expint)
uppergamma, sinh, atan, Ne, expint, Integral)

from sympy import beta as beta_fn
from sympy import cos, sin, tan, atan, exp, besseli, besselj, besselk
Expand Down Expand Up @@ -1563,58 +1563,78 @@ def GammaInverse(name, a, b):
return rv(name, GammaInverseDistribution, (a, b))

#-------------------------------------------------------------------------------
# Gumbel distribution --------------------------------------------------------
# Gumbel distribution (Maximum and Minimum) --------------------------------------------------------


class GumbelDistribution(SingleContinuousDistribution):
_argnames = ('beta', 'mu')
_argnames = ('beta', 'mu', 'minimum')

set = Interval(-oo, oo)

@staticmethod
def check(beta, mu):
def check(beta, mu, minimum):
_value_check(beta > 0, "Scale parameter beta must be positive.")

def pdf(self, x):
beta, mu = self.beta, self.mu
z = (x - mu)/beta
return (1/beta)*exp(-(z + exp(-z)))
f_max = (1/beta)*exp(-z - exp(-z))
f_min = (1/beta)*exp(z - exp(z))
return Piecewise((f_min, self.minimum), (f_max, not self.minimum))

def _cdf(self, x):
beta, mu = self.beta, self.mu
return exp(-exp((mu - x)/beta))
z = (x - mu)/beta
F_max = exp(-exp(-z))
F_min = 1 - exp(-exp(z))
return Piecewise((F_min, self.minimum), (F_max, not self.minimum))

def _characteristic_function(self, t):
return gamma(1 - I*self.beta*t) * exp(I*self.mu*t)
cf_max = gamma(1 - I*self.beta*t) * exp(I*self.mu*t)
cf_min = gamma(1 + I*self.beta*t) * exp(I*self.mu*t)
return Piecewise((cf_min, self.minimum), (cf_max, not self.minimum))

def _moment_generating_function(self, t):
return gamma(1 - self.beta*t) * exp(self.mu*t)
mgf_max = gamma(1 - self.beta*t) * exp(self.mu*t)
mgf_min = gamma(1 + self.beta*t) * exp(self.mu*t)
return Piecewise((mgf_min, self.minimum), (mgf_max, not self.minimum))

def Gumbel(name, beta, mu):
def Gumbel(name, beta, mu, minimum=False):
r"""
Create a Continuous Random Variable with Gumbel distribution.

The density of the Gumbel distribution is given by

For Maximum

.. math::
f(x) := \dfrac{1}{\beta} \exp \left( -\dfrac{x-\mu}{\beta}
- \exp \left( -\dfrac{x - \mu}{\beta} \right) \right)

with :math:`x \in [ - \infty, \infty ]`.

For Minimum

.. math::
f(x) := \frac{e^{- e^{\frac{- \mu + x}{\beta}} + \frac{- \mu + x}{\beta}}}{\beta}

with :math:`x \in [ - \infty, \infty ]`.

Parameters
==========

mu: Real number, 'mu' is a location
beta: Real number, 'beta > 0' is a scale
mu : Real number, 'mu' is a location
beta : Real number, 'beta > 0' is a scale
minimum : Boolean, by default, False, set to True for enabling minimum distribution

Returns
==========
=======

A RandomSymbol

Examples
==========
========

>>> from sympy.stats import Gumbel, density, E, variance, cdf
>>> from sympy import Symbol, simplify, pprint
>>> x = Symbol("x")
Expand All @@ -1624,16 +1644,18 @@ def Gumbel(name, beta, mu):
>>> density(X)(x)
exp(-exp(-(-mu + x)/beta) - (-mu + x)/beta)/beta
>>> cdf(X)(x)
exp(-exp((mu - x)/beta))
exp(-exp(-(-mu + x)/beta))

References
==========

.. [1] http://mathworld.wolfram.com/GumbelDistribution.html
.. [2] https://en.wikipedia.org/wiki/Gumbel_distribution
.. [3] http://www.mathwave.com/help/easyfit/html/analyses/distributions/gumbel_max.html
.. [4] http://www.mathwave.com/help/easyfit/html/analyses/distributions/gumbel_min.html

"""
return rv(name, GumbelDistribution, (beta, mu))
return rv(name, GumbelDistribution, (beta, mu, minimum))

#-------------------------------------------------------------------------------
# Gompertz distribution --------------------------------------------------------
Expand Down
10 changes: 8 additions & 2 deletions sympy/stats/tests/test_continuous_rv.py
Expand Up @@ -631,9 +631,15 @@ def test_gumbel():
beta = Symbol("beta", positive=True)
mu = Symbol("mu")
x = Symbol("x")
y = Symbol("y")
X = Gumbel("x", beta, mu)
assert str(density(X)(x)) == 'exp(-exp(-(-mu + x)/beta) - (-mu + x)/beta)/beta'
assert cdf(X)(x) == exp(-exp((mu - x)/beta))
Y = Gumbel("y", beta, mu, minimum=True)
assert density(X)(x).expand() == \
exp(mu/beta)*exp(-x/beta)*exp(-exp(mu/beta)*exp(-x/beta))/beta
assert density(Y)(y).expand() == \
exp(-mu/beta)*exp(y/beta)*exp(-exp(-mu/beta)*exp(y/beta))/beta
assert cdf(X)(x).expand() == \
exp(-exp(mu/beta)*exp(-x/beta))


def test_kumaraswamy():
Expand Down