Skip to content
This repository

WIP: More distributions #1413

Merged
merged 17 commits into from over 1 year ago

7 participants

raoulb Don't Add Me To Your Organization a.k.a The Travis Bot Stefan Krastanov Ondřej Čertík Aaron Meurer Julien Rioux Matthew Rocklin
raoulb
Collaborator
raoulb commented

Even more continuous distribution functions

Don't Add Me To Your Organization a.k.a The Travis Bot

This pull request fails (merged 26c4adfc into 9ff4edf).

Stefan Krastanov
Collaborator

SymPy Bot Summary: :red_circle: There were test failures.

@raoulb: Please fix the test failures.

Test command: setup.py test
master hash: 9ff4edf
branch hash: 26c4adf

Interpreter 1: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY_-0fDA

Interpreter 2: :red_circle: There were test failures.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYxq8fDA

Interpreter 3: :red_circle: There were test failures.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYz9YfDA

Build HTML Docs: :red_circle: There were test failures.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYvYgfDA

Automatic review by SymPy Bot.

Ondřej Čertík
Owner

It needs regular tests before we can merge it.

raoulb
Collaborator

I plan to add some more distribution in this PR.
Let's just keep it open for the moment.

Don't Add Me To Your Organization a.k.a The Travis Bot

This pull request fails (merged 8691d70b into a184841).

Don't Add Me To Your Organization a.k.a The Travis Bot

This pull request passes (merged e49899e8 into a184841).

raoulb
Collaborator

Do we want this to get into SymPy 0.7.2?
If yes, so lets merge it now.
If no, I will continue adding distributions on top of this PR some day.

Don't Add Me To Your Organization a.k.a The Travis Bot

This pull request passes (merged 7508eec6 into a184841).

Aaron Meurer
Owner

I already created the 0.7.2 branch and would like to keep it for bug fixes and API changes only, so let's shoot for 0.7.3 for this.

Even so, if the work here is ok as is, this can be merged into master and you can open new PRs for further distributions.

Julien Rioux
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: All tests have passed.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYhYYlDA

Interpreter: /usr/bin/python (2.7.0-final-0)
Architecture: Linux (32-bit)
Cache: yes
Test command: setup.py test
master hash: 161c963
branch hash: 7508eec

Automatic review by SymPy Bot.

Julien Rioux
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: All tests have passed.

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYkOIjDA

Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 161c963
branch hash: 7508eec

Automatic review by SymPy Bot.

Julien Rioux
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: Passed after merging raoulb/more_distributions (7508eec) into master (c49dca1).
:eight_spoked_asterisk:Python 2.7.2-final-0: pass
:eight_spoked_asterisk:Python 3.2.1-final-0: pass
:eight_spoked_asterisk:Sphinx 1.1.3: pass
Docs build command: make clean && make html-errors && make clean && make latex && cd _build/latex && xelatex -interaction=nonstopmode sympy-*.tex

Julien Rioux
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: Passed after merging raoulb/more_distributions (a24f147) into master (93f6ce3).
:eight_spoked_asterisk:Python 2.7.2-final-0: pass
:eight_spoked_asterisk:Python 3.2.1-final-0: pass
:eight_spoked_asterisk:Sphinx 1.1.3: pass
Docs build command: make clean && make html-errors && make clean && make latex && cd _build/latex && xelatex sympy-*.tex

Julien Rioux
Collaborator

SymPy Bot Summary: :exclamation: There were merge conflicts (could not merge raoulb/more_distributions (a24f147) into master (57e94e4)); could not test the branch.
@raoulb: Please rebase or merge your branch with master. See the report for a list of the merge conflicts.

raoulb
Collaborator

Is this ready for merge?

sympy/stats/tests/test_continuous_rv.py
... ... @@ -1,14 +1,18 @@
1 1 from sympy.stats import (P, E, where, density, variance, covariance, skewness,
2   - given, pspace, cdf, ContinuousRV, sample,
3   - Arcsin, Benini, Beta, BetaPrime, Cauchy, Chi, Dagum,
4   - Exponential, Gamma, Laplace, Logistic, LogNormal,
5   - Maxwell, Nakagami, Normal, Pareto, Rayleigh,
6   - StudentT, Triangular, Uniform, UniformSum, Weibull,
7   - WignerSemicircle)
  2 + given, pspace, cdf, ContinuousRV, sample
1
Matthew Rocklin Collaborator
mrocklin added a note

Missing a comma here.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Matthew Rocklin
Collaborator

I'm playing with sympy.stats again. It'd be nice to get this in before I make substantial changes. I've commented next to a trivial syntax error in the test file.

Matthew Rocklin
Collaborator

Also, I'm curious if there was a reason why you haven't implemented ChiSquared? You have Chi and ChiNonCentral but my understanding is that ChiSquared is used more often than either.

Matthew Rocklin
Collaborator

I ask just because I'm about to do it and was wondering if there was a reason why it was skipped. (My knowledge of distributions is very small).

raoulb
Collaborator

Also, I'm curious if there was a reason why you haven't implemented ChiSquared? You have Chi and ChiNonCentral but my understanding is that ChiSquared is used more often than either.

I think the reason was we can get if from \chi easily: let X ~ \chi^2_k then \sqrt{X} ~ \chi_k.
But maybe you are right and we should add it too, at least for the reason that users do not miss it.

raoulb
Collaborator

I'm playing with sympy.stats again. It'd be nice to get this in before I make substantial changes.

Yes, let's merge this one.
What are your new plans for the stats module?

BTW: The pattern/rule based RV-simplifier is a great project! But I doubt that I find time working on it.

Matthew Rocklin
Collaborator

The simplifier is the main point. Sympy.stats will require some general refactoring before it's ready though.

Sorry to hear you're busy. You're quite a force. I still hold out hope that I'll be able to attract you to write down rewrite rules once they're easy to do.

I'll merge in a bit.

sympy/stats/crv_types.py
((88 lines not shown))
  598 + ==========
  599 +
  600 + k : Integer, `k` > 0 the number of degrees of freedom
  601 +
  602 + Returns
  603 + =======
  604 +
  605 + A RandomSymbol.
  606 +
  607 + Examples
  608 + ========
  609 +
  610 + >>> from sympy.stats import ChiSquared, density, E, std
  611 + >>> from sympy import Symbol, simplify, combsimp
  612 +
  613 + >>> k = Symbol("k", integer=True)
1
Matthew Rocklin Collaborator
mrocklin added a note

Looks like we might need positive=True as well

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/stats/crv_types.py
((99 lines not shown))
  609 +
  610 + >>> from sympy.stats import ChiSquared, density, E, std
  611 + >>> from sympy import Symbol, simplify, combsimp
  612 +
  613 + >>> k = Symbol("k", integer=True)
  614 +
  615 + >>> X = ChiSquared("x", k)
  616 +
  617 + >>> density(X)
  618 + Lambda(_x, 2**(-k/2)*_x**(k/2 - 1)*exp(-_x/2)/gamma(k/2))
  619 +
  620 + >>> combsimp(E(X))
  621 + k
  622 +
  623 + >>> simplify(expand_func(variance(X)))
  624 + >>> 2*k
1
Matthew Rocklin Collaborator
mrocklin added a note

remove >>>

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Matthew Rocklin mrocklin commented on the diff
sympy/stats/crv_types.py
((98 lines not shown))
  608 + ========
  609 +
  610 + >>> from sympy.stats import ChiSquared, density, E, std
  611 + >>> from sympy import Symbol, simplify, combsimp
  612 +
  613 + >>> k = Symbol("k", integer=True)
  614 +
  615 + >>> X = ChiSquared("x", k)
  616 +
  617 + >>> density(X)
  618 + Lambda(_x, 2**(-k/2)*_x**(k/2 - 1)*exp(-_x/2)/gamma(k/2))
  619 +
  620 + >>> combsimp(E(X))
  621 + k
  622 +
  623 + >>> simplify(expand_func(variance(X)))
2
Matthew Rocklin Collaborator
mrocklin added a note

expand_func needs to be imported

raoulb Collaborator
raoulb added a note

Oh, right. I had to add it when I tried this in ipython first and then forgot to put it into the docstring as well.
It turned out I forgot to run the doctests, there are some more missing imports and typos.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
raoulb
Collaborator

SymPy Bot Summary: :red_circle: Failed after merging raoulb/more_distributions (4f3ea61) into master (623aead).
@raoulb: Please fix the test failures.
:red_circle:Python 2.7.3-final-0: fail

raoulb
Collaborator

Some more unrelated bot failures

raoulb
Collaborator

Let's wait for travis and then merge

Matthew Rocklin mrocklin merged commit 143860f into from
Matthew Rocklin
Collaborator

Merged.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
11 doc/src/modules/stats.rst
Source Rendered
@@ -25,9 +25,17 @@ Continuous Types
25 25 .. autofunction:: BetaPrime
26 26 .. autofunction:: Cauchy
27 27 .. autofunction:: Chi
  28 +.. autofunction:: ChiNoncentral
  29 +.. autofunction:: ChiSquared
28 30 .. autofunction:: Dagum
  31 +.. autofunction:: Erlang
29 32 .. autofunction:: Exponential
  33 +.. autofunction:: FDistribution
  34 +.. autofunction:: FisherZ
  35 +.. autofunction:: Frechet
30 36 .. autofunction:: Gamma
  37 +.. autofunction:: GammaInverse
  38 +.. autofunction:: Kumaraswamy
31 39 .. autofunction:: Laplace
32 40 .. autofunction:: Logistic
33 41 .. autofunction:: LogNormal
@@ -35,11 +43,14 @@ Continuous Types
35 43 .. autofunction:: Nakagami
36 44 .. autofunction:: Normal
37 45 .. autofunction:: Pareto
  46 +.. autofunction:: QuadraticU
  47 +.. autofunction:: RaisedCosine
38 48 .. autofunction:: Rayleigh
39 49 .. autofunction:: StudentT
40 50 .. autofunction:: Triangular
41 51 .. autofunction:: Uniform
42 52 .. autofunction:: UniformSum
  53 +.. autofunction:: VonMises
43 54 .. autofunction:: Weibull
44 55 .. autofunction:: WignerSemicircle
45 56 .. autofunction:: ContinuousRV
47 sympy/core/tests/test_args.py
@@ -696,6 +696,16 @@ def test_sympy__stats__crv_types__ChiPSpace():
696 696 assert _test_args(ChiPSpace('X', 1))
697 697
698 698
  699 +def test_sympy__stats__crv_types__ChiNoncentralPSpace():
  700 + from sympy.stats.crv_types import ChiNoncentralPSpace
  701 + assert _test_args(ChiNoncentralPSpace('X', 1,1))
  702 +
  703 +
  704 +def test_sympy__stats__crv_types__ChiSquaredPSpace():
  705 + from sympy.stats.crv_types import ChiSquaredPSpace
  706 + assert _test_args(ChiSquaredPSpace('X', 1))
  707 +
  708 +
699 709 def test_sympy__stats__crv_types__DagumPSpace():
700 710 from sympy.stats.crv_types import DagumPSpace
701 711 assert _test_args(DagumPSpace('X', 1, 1, 1))
@@ -706,11 +716,35 @@ def test_sympy__stats__crv_types__ExponentialPSpace():
706 716 assert _test_args(ExponentialPSpace('X', 1))
707 717
708 718
  719 +def test_sympy__stats__crv_types__FDistributionPSpace():
  720 + from sympy.stats.crv_types import FDistributionPSpace
  721 + assert _test_args(FDistributionPSpace('X', 1, 1))
  722 +
  723 +
  724 +def test_sympy__stats__crv_types__FisherZPSpace():
  725 + from sympy.stats.crv_types import FisherZPSpace
  726 + assert _test_args(FisherZPSpace('X', 1, 1))
  727 +
  728 +
  729 +def test_sympy__stats__crv_types__FrechetPSpace():
  730 + from sympy.stats.crv_types import FrechetPSpace
  731 + assert _test_args(FrechetPSpace('X', 1, 1, 1))
  732 +
  733 +
  734 +def test_sympy__stats__crv_types__GammaInversePSpace():
  735 + from sympy.stats.crv_types import GammaInversePSpace
  736 + assert _test_args(GammaInversePSpace('X', 1, 1))
  737 +
  738 +
709 739 def test_sympy__stats__crv_types__GammaPSpace():
710 740 from sympy.stats.crv_types import GammaPSpace
711 741 assert _test_args(GammaPSpace('X', 1, 1))
712 742
713 743
  744 +def test_sympy__stats__crv_types__KumaraswamyPSpace():
  745 + from sympy.stats.crv_types import KumaraswamyPSpace
  746 + assert _test_args(KumaraswamyPSpace('X', 1, 1))
  747 +
714 748 def test_sympy__stats__crv_types__LaplacePSpace():
715 749 from sympy.stats.crv_types import LaplacePSpace
716 750 assert _test_args(LaplacePSpace('X', 0, 1))
@@ -746,6 +780,14 @@ def test_sympy__stats__crv_types__ParetoPSpace():
746 780 assert _test_args(ParetoPSpace('X', 1, 1))
747 781
748 782
  783 +def test_sympy__stats__crv_types__QuadraticUPSpace():
  784 + from sympy.stats.crv_types import QuadraticUPSpace
  785 + assert _test_args(QuadraticUPSpace('X', 1, 2))
  786 +
  787 +def test_sympy__stats__crv_types__RaisedCosinePSpace():
  788 + from sympy.stats.crv_types import RaisedCosinePSpace
  789 + assert _test_args(RaisedCosinePSpace('X', 1, 1))
  790 +
749 791 def test_sympy__stats__crv_types__RayleighPSpace():
750 792 from sympy.stats.crv_types import RayleighPSpace
751 793 assert _test_args(RayleighPSpace('X', 1))
@@ -771,6 +813,11 @@ def test_sympy__stats__crv_types__UniformSumPSpace():
771 813 assert _test_args(UniformSumPSpace('X', 1))
772 814
773 815
  816 +def test_sympy__stats__crv_types__VonMisesPSpace():
  817 + from sympy.stats.crv_types import VonMisesPSpace
  818 + assert _test_args(VonMisesPSpace('X', 1, 1))
  819 +
  820 +
774 821 def test_sympy__stats__crv_types__WeibullPSpace():
775 822 from sympy.stats.crv_types import WeibullPSpace
776 823 assert _test_args(WeibullPSpace('X', 1, 1))
11 sympy/stats/__init__.py
@@ -53,9 +53,12 @@
53 53
54 54 import crv_types
55 55 from crv_types import (
56   - Arcsin, Benini, Beta, BetaPrime, Cauchy, Chi, ContinuousRV, Dagum,
57   - Exponential, Gamma, Laplace, Logistic, LogNormal, Maxwell, Nakagami,
58   - Normal, Pareto, Rayleigh, StudentT, Triangular, Uniform, UniformSum,
59   - Weibull, WignerSemicircle,
  56 + ContinuousRV,
  57 + Arcsin, Benini, Beta, BetaPrime, Cauchy, Chi, ChiNoncentral, ChiSquared,
  58 + Dagum, Erlang, Exponential, FDistribution, FisherZ, Frechet, Gamma,
  59 + GammaInverse, Kumaraswamy, Laplace, Logistic, LogNormal, Maxwell,
  60 + Nakagami, Normal, Pareto, QuadraticU, RaisedCosine, Rayleigh,
  61 + StudentT, Triangular, Uniform, UniformSum, VonMises, Weibull,
  62 + WignerSemicircle
60 63 )
61 64 __all__.extend(crv_types.__all__)
752 sympy/stats/crv_types.py
@@ -9,9 +9,17 @@
9 9 BetaPrime
10 10 Cauchy
11 11 Chi
  12 +ChiNoncentral
  13 +ChiSquared
12 14 Dagum
  15 +Erlang
13 16 Exponential
  17 +FDistribution
  18 +FisherZ
  19 +Frechet
14 20 Gamma
  21 +GammaInverse
  22 +Kumaraswamy
15 23 Laplace
16 24 Logistic
17 25 LogNormal
@@ -19,19 +27,23 @@
19 27 Nakagami
20 28 Normal
21 29 Pareto
  30 +QuadraticU
  31 +RaisedCosine
22 32 Rayleigh
23 33 StudentT
24 34 Triangular
25 35 Uniform
26 36 UniformSum
  37 +VonMises
27 38 Weibull
28 39 WignerSemicircle
29 40 """
30 41
31 42 from sympy import (exp, log, sqrt, pi, S, Dummy, Interval, S, sympify, gamma,
32 43 Piecewise, And, Eq, binomial, factorial, Sum, floor, Abs,
33   - Symbol, log)
  44 + Symbol, log, besseli)
34 45 from sympy import beta as beta_fn
  46 +from sympy import cos, exp, besseli
35 47 from crv import SingleContinuousPSpace
36 48 from sympy.core.decorators import _sympifyit
37 49 import random
@@ -45,9 +57,17 @@
45 57 'BetaPrime',
46 58 'Cauchy',
47 59 'Chi',
  60 +'ChiNoncentral',
  61 +'ChiSquared',
48 62 'Dagum',
  63 +'Erlang',
49 64 'Exponential',
  65 +'FDistribution',
  66 +'FisherZ',
  67 +'Frechet',
50 68 'Gamma',
  69 +'GammaInverse',
  70 +'Kumaraswamy',
51 71 'Laplace',
52 72 'Logistic',
53 73 'LogNormal',
@@ -55,11 +75,14 @@
55 75 'Nakagami',
56 76 'Normal',
57 77 'Pareto',
  78 +'QuadraticU',
  79 +'RaisedCosine',
58 80 'Rayleigh',
59 81 'StudentT',
60 82 'Triangular',
61 83 'Uniform',
62 84 'UniformSum',
  85 +'VonMises',
63 86 'Weibull',
64 87 'WignerSemicircle'
65 88 ]
@@ -488,6 +511,128 @@ def Chi(name, k):
488 511 return ChiPSpace(name, k).value
489 512
490 513 #-------------------------------------------------------------------------------
  514 +# Non-central Chi distribution -------------------------------------------------
  515 +
  516 +
  517 +class ChiNoncentralPSpace(SingleContinuousPSpace):
  518 + def __new__(cls, name, k, l):
  519 + k = sympify(k)
  520 + l = sympify(l)
  521 + x = Symbol(name)
  522 + pdf = exp(-(x**2+l**2)/2)*x**k*l / (l*x)**(k/2) * besseli(k/2-1, l*x)
  523 + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set = Interval(0, oo))
  524 + return obj
  525 +
  526 +
  527 +def ChiNoncentral(name, k, l):
  528 + r"""
  529 + Create a continuous random variable with a non-central Chi distribution.
  530 +
  531 + The density of the non-central Chi distribution is given by
  532 +
  533 + .. math::
  534 + f(x) := \frac{e^{-(x^2+\lambda^2)/2} x^k\lambda}
  535 + {(\lambda x)^{k/2}} I_{k/2-1}(\lambda x)
  536 +
  537 + with :math:`x \geq 0`.
  538 +
  539 + Parameters
  540 + ==========
  541 +
  542 + k : `k` > 0 the number of degrees of freedom
  543 + l : Shift parameter
  544 +
  545 + Returns
  546 + =======
  547 +
  548 + A RandomSymbol.
  549 +
  550 + Examples
  551 + ========
  552 +
  553 + >>> from sympy.stats import ChiNoncentral, density, E, std
  554 + >>> from sympy import Symbol, simplify
  555 +
  556 + >>> k = Symbol("k", integer=True)
  557 + >>> l = Symbol("l")
  558 +
  559 + >>> X = ChiNoncentral("x", k, l)
  560 +
  561 + >>> density(X)
  562 + Lambda(_x, _x**k*l*(_x*l)**(-k/2)*exp(-_x**2/2 - l**2/2)*besseli(k/2 - 1, _x*l))
  563 +
  564 + References
  565 + ==========
  566 +
  567 + [1] http://en.wikipedia.org/wiki/Noncentral_chi_distribution
  568 + """
  569 +
  570 + return ChiNoncentralPSpace(name, k, l).value
  571 +
  572 +#-------------------------------------------------------------------------------
  573 +# Chi squared distribution -----------------------------------------------------
  574 +
  575 +
  576 +class ChiSquaredPSpace(SingleContinuousPSpace):
  577 + def __new__(cls, name, k):
  578 + k = sympify(k)
  579 + x = Symbol(name)
  580 + pdf = 1/(2**(k/2)*gamma(k/2))*x**(k/2 - 1)*exp(-x/2)
  581 + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, oo))
  582 + return obj
  583 +
  584 +
  585 +def ChiSquared(name, k):
  586 + r"""
  587 + Create a continuous random variable with a Chi-squared distribution.
  588 +
  589 + The density of the Chi-squared distribution is given by
  590 +
  591 + .. math::
  592 + f(x) := \frac{1}{2^{\frac{k}{2}}\Gamma\left(\frac{k}{2}\right)}
  593 + x^{\frac{k}{2}-1} e^{-\frac{x}{2}}
  594 +
  595 + with :math:`x \geq 0`.
  596 +
  597 + Parameters
  598 + ==========
  599 +
  600 + k : Integer, `k` > 0 the number of degrees of freedom
  601 +
  602 + Returns
  603 + =======
  604 +
  605 + A RandomSymbol.
  606 +
  607 + Examples
  608 + ========
  609 +
  610 + >>> from sympy.stats import ChiSquared, density, E, variance
  611 + >>> from sympy import Symbol, simplify, combsimp, expand_func
  612 +
  613 + >>> k = Symbol("k", integer=True, positive=True)
  614 +
  615 + >>> X = ChiSquared("x", k)
  616 +
  617 + >>> density(X)
  618 + Lambda(_x, 2**(-k/2)*_x**(k/2 - 1)*exp(-_x/2)/gamma(k/2))
  619 +
  620 + >>> combsimp(E(X))
  621 + k
  622 +
  623 + >>> simplify(expand_func(variance(X)))
  624 + 2*k
  625 +
  626 + References
  627 + ==========
  628 +
  629 + [1] http://en.wikipedia.org/wiki/Chi_squared_distribution
  630 + [2] http://mathworld.wolfram.com/Chi-SquaredDistribution.html
  631 + """
  632 +
  633 + return ChiSquaredPSpace(name, k).value
  634 +
  635 +#-------------------------------------------------------------------------------
491 636 # Dagum distribution -----------------------------------------------------------
492 637
493 638
@@ -548,6 +693,72 @@ def Dagum(name, p, a, b):
548 693 return DagumPSpace(name, p, a, b).value
549 694
550 695 #-------------------------------------------------------------------------------
  696 +# Erlang distribution ----------------------------------------------------------
  697 +
  698 +def Erlang(name, k, l):
  699 + r"""
  700 + Create a continuous random variable with an Erlang distribution.
  701 +
  702 + The density of the Erlang distribution is given by
  703 +
  704 + .. math::
  705 + f(x) := \frac{\lambda^k x^{k-1} e^{-\lambda x}}{(k-1)!}
  706 +
  707 + with :math:`x \in [0,\infty]`.
  708 +
  709 + Parameters
  710 + ==========
  711 +
  712 + k : Integer
  713 + l : Real number, :math:`\lambda` > 0 the rate
  714 +
  715 + Returns
  716 + =======
  717 +
  718 + A RandomSymbol.
  719 +
  720 + Examples
  721 + ========
  722 +
  723 + >>> from sympy.stats import Erlang, density, cdf, E, variance
  724 + >>> from sympy import Symbol, simplify, pprint
  725 +
  726 + >>> k = Symbol("k", integer=True, positive=True)
  727 + >>> l = Symbol("l", positive=True)
  728 +
  729 + >>> X = Erlang("x", k, l)
  730 +
  731 + >>> D = density(X)
  732 + >>> pprint(D, use_unicode=False)
  733 + / k - 1 k -x*l\
  734 + | x *l *e |
  735 + Lambda|x, ---------------|
  736 + \ gamma(k) /
  737 +
  738 + >>> C = cdf(X, meijerg=True)
  739 + >>> pprint(C, use_unicode=False)
  740 + / / k*lowergamma(k, 0) k*lowergamma(k, z*l) \
  741 + Lambda|z, |- ------------------ + -------------------- for z >= 0|
  742 + | < gamma(k + 1) gamma(k + 1) |
  743 + | | |
  744 + \ \ 0 otherwise /
  745 +
  746 + >>> simplify(E(X))
  747 + k/l
  748 +
  749 + >>> simplify(variance(X))
  750 + (gamma(k)*gamma(k + 2) - gamma(k + 1)**2)/(l**2*gamma(k)**2)
  751 +
  752 + References
  753 + ==========
  754 +
  755 + [1] http://en.wikipedia.org/wiki/Erlang_distribution
  756 + [2] http://mathworld.wolfram.com/ErlangDistribution.html
  757 + """
  758 +
  759 + return GammaPSpace(name, k, 1/l).value
  760 +
  761 +#-------------------------------------------------------------------------------
551 762 # Exponential distribution -----------------------------------------------------
552 763
553 764
@@ -636,6 +847,202 @@ def Exponential(name, rate):
636 847 return ExponentialPSpace(name, rate).value
637 848
638 849 #-------------------------------------------------------------------------------
  850 +# F distribution ---------------------------------------------------------------
  851 +
  852 +class FDistributionPSpace(SingleContinuousPSpace):
  853 + def __new__(cls, name, d1, d2):
  854 + d1 = sympify(d1)
  855 + d2 = sympify(d2)
  856 + x = Symbol(name)
  857 + pdf = (sqrt((d1*x)**d1*d2**d2 / (d1*x+d2)**(d1+d2))
  858 + / (x * beta_fn(d1/2, d2/2)))
  859 + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, oo))
  860 + return obj
  861 +
  862 +def FDistribution(name, d1, d2):
  863 + r"""
  864 + Create a continuous random variable with a F distribution.
  865 +
  866 + The density of the F distribution is given by
  867 +
  868 + .. math::
  869 + f(x) := \frac{\sqrt{\frac{(d_1 x)^{d_1} d_2^{d_2}}
  870 + {(d_1 x + d_2)^{d_1 + d_2}}}}
  871 + {x \mathrm{B} \left(\frac{d_1}{2}, \frac{d_2}{2}\right)}
  872 +
  873 + with :math:`x > 0`.
  874 +
  875 + Parameters
  876 + ==========
  877 +
  878 + d1 : `d1` > 0 a parameter
  879 + d2 : `d2` > 0 a parameter
  880 +
  881 + Returns
  882 + =======
  883 +
  884 + A RandomSymbol.
  885 +
  886 + Examples
  887 + ========
  888 +
  889 + >>> from sympy.stats import FDistribution, density
  890 + >>> from sympy import Symbol, simplify, pprint
  891 +
  892 + >>> d1 = Symbol("d1", positive=True)
  893 + >>> d2 = Symbol("d2", positive=True)
  894 +
  895 + >>> X = FDistribution("x", d1, d2)
  896 +
  897 + >>> D = density(X)
  898 + >>> pprint(D, use_unicode=False)
  899 + / d2 \
  900 + | -- ______________________________ |
  901 + | 2 / d1 -d1 - d2 /d1 d2\|
  902 + | d2 *\/ (x*d1) *(x*d1 + d2) *gamma|-- + --||
  903 + | \2 2 /|
  904 + Lambda|x, -----------------------------------------------------|
  905 + | /d1\ /d2\ |
  906 + | x*gamma|--|*gamma|--| |
  907 + \ \2 / \2 / /
  908 +
  909 + References
  910 + ==========
  911 +
  912 + [1] http://en.wikipedia.org/wiki/F-distribution
  913 + [2] http://mathworld.wolfram.com/F-Distribution.html
  914 + """
  915 +
  916 + return FDistributionPSpace(name, d1, d2).value
  917 +
  918 +#-------------------------------------------------------------------------------
  919 +# Fisher Z distribution --------------------------------------------------------
  920 +
  921 +class FisherZPSpace(SingleContinuousPSpace):
  922 + def __new__(cls, name, d1, d2):
  923 + d1 = sympify(d1)
  924 + d2 = sympify(d2)
  925 + x = Symbol(name)
  926 + pdf = (2*d1**(d1/2)*d2**(d2/2) / beta_fn(d1/2, d2/2) *
  927 + exp(d1*x) / (d1*exp(2*x)+d2)**((d1+d2)/2))
  928 + obj = SingleContinuousPSpace.__new__(cls, x, pdf)
  929 + return obj
  930 +
  931 +def FisherZ(name, d1, d2):
  932 + r"""
  933 + Create a Continuous Random Variable with an Fisher's Z distribution.
  934 +
  935 + The density of the Fisher's Z distribution is given by
  936 +
  937 + .. math::
  938 + f(x) := \frac{2d_1^{d_1/2} d_2^{d_2/2}} {\mathrm{B}(d_1/2, d_2/2)}
  939 + \frac{e^{d_1z}}{\left(d_1e^{2z}+d_2\right)^{\left(d_1+d_2\right)/2}}
  940 +
  941 + Parameters
  942 + ==========
  943 +
  944 + d1 : `d1` > 0, degree of freedom
  945 + d2 : `d2` > 0, degree of freedom
  946 +
  947 + Returns
  948 + =======
  949 +
  950 + A RandomSymbol.
  951 +
  952 + Examples
  953 + ========
  954 +
  955 + >>> from sympy.stats import FisherZ, density
  956 + >>> from sympy import Symbol, simplify, pprint
  957 +
  958 + >>> d1 = Symbol("d1", positive=True)
  959 + >>> d2 = Symbol("d2", positive=True)
  960 +
  961 + >>> X = FisherZ("x", d1, d2)
  962 +
  963 + >>> D = density(X)
  964 + >>> pprint(D, use_unicode=False)
  965 + / d1 d2 \
  966 + | d1 d2 - -- - -- |
  967 + | -- -- 2 2 |
  968 + | 2 2 / 2*x \ x*d1 /d1 d2\|
  969 + | 2*d1 *d2 *\d1*e + d2/ *e *gamma|-- + --||
  970 + | \2 2 /|
  971 + Lambda|x, --------------------------------------------------------|
  972 + | /d1\ /d2\ |
  973 + | gamma|--|*gamma|--| |
  974 + \ \2 / \2 / /
  975 +
  976 + References
  977 + ==========
  978 +
  979 + [1] http://en.wikipedia.org/wiki/Fisher%27s_z-distribution
  980 + [2] http://mathworld.wolfram.com/Fishersz-Distribution.html
  981 + """
  982 +
  983 + return FisherZPSpace(name, d1, d2).value
  984 +
  985 +#-------------------------------------------------------------------------------
  986 +# Frechet distribution ---------------------------------------------------------
  987 +
  988 +class FrechetPSpace(SingleContinuousPSpace):
  989 + def __new__(cls, name, a, s=0, m=0):
  990 + a = sympify(a)
  991 + s = sympify(s)
  992 + m = sympify(m)
  993 + x = Symbol(name)
  994 + pdf = a/s * ((x-m)/s)**(-1-a) * exp(-((x-m)/s)-a)
  995 + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set = Interval(m, oo))
  996 + return obj
  997 +
  998 +def Frechet(name, a, s=1, m=0):
  999 + r"""
  1000 + Create a continuous random variable with a Frechet distribution.
  1001 +
  1002 + The density of the Frechet distribution is given by
  1003 +
  1004 + .. math::
  1005 + f(x) := \frac{\alpha}{s} \left(\frac{x-m}{s}\right)^{-1-\alpha}
  1006 + e^{-(\frac{x-m}{s})^{-\alpha}}
  1007 +
  1008 + with :math:`x \geq m`.
  1009 +
  1010 + Parameters
  1011 + ==========
  1012 +
  1013 + a : Real number, :math:`a \in \left(0, \infty\right)` the shape
  1014 + s : Real number, :math:`s \in \left(0, \infty\right)` the scale
  1015 + m : Real number, :math:`m \in \left(-\infty, \infty\right)` the minimum
  1016 +
  1017 + Returns
  1018 + =======
  1019 +
  1020 + A RandomSymbol.
  1021 +
  1022 + Examples
  1023 + ========
  1024 +
  1025 + >>> from sympy.stats import Frechet, density, E, std
  1026 + >>> from sympy import Symbol, simplify
  1027 +
  1028 + >>> a = Symbol("a", positive=True)
  1029 + >>> s = Symbol("s", positive=True)
  1030 + >>> m = Symbol("m", real=True)
  1031 +
  1032 + >>> X = Frechet("x", a, s, m)
  1033 +
  1034 + >>> density(X)
  1035 + Lambda(_x, a*((_x - m)/s)**(-a - 1)*exp(-a - (_x - m)/s)/s)
  1036 +
  1037 + References
  1038 + ==========
  1039 +
  1040 + [1] http://en.wikipedia.org/wiki/Fr%C3%A9chet_distribution
  1041 + """
  1042 +
  1043 + return FrechetPSpace(name, a, s, m).value
  1044 +
  1045 +#-------------------------------------------------------------------------------
639 1046 # Gamma distribution -----------------------------------------------------------
640 1047
641 1048
@@ -731,6 +1138,141 @@ def Gamma(name, k, theta):
731 1138 return GammaPSpace(name, k, theta).value
732 1139
733 1140 #-------------------------------------------------------------------------------
  1141 +# Inverse Gamma distribution ---------------------------------------------------
  1142 +
  1143 +class GammaInversePSpace(SingleContinuousPSpace):
  1144 + def __new__(cls, name, a, b):
  1145 + a = sympify(a)
  1146 + b = sympify(b)
  1147 + _value_check(a > 0, "alpha must be positive")
  1148 + _value_check(b > 0, "beta must be positive")
  1149 + x = Symbol(name)
  1150 + pdf = b**a/gamma(a) * x**(-a-1) * exp(-b/x)
  1151 + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, oo))
  1152 + obj.a = a
  1153 + obj.b = b
  1154 + return obj
  1155 +
  1156 +def GammaInverse(name, a, b):
  1157 + r"""
  1158 + Create a continuous random variable with an inverse Gamma distribution.
  1159 +
  1160 + The density of the inverse Gamma distribution is given by
  1161 +
  1162 + .. math::
  1163 + f(x) := \frac{\beta^\alpha}{\Gamma(\alpha)} x^{-\alpha - 1}
  1164 + \exp\left(\frac{-\beta}{x}\right)
  1165 +
  1166 + with :math:`x > 0`.
  1167 +
  1168 + Parameters
  1169 + ==========
  1170 +
  1171 + a : Real number, `a` > 0 a shape
  1172 + b : Real number, `b` > 0 a scale
  1173 +
  1174 + Returns
  1175 + =======
  1176 +
  1177 + A RandomSymbol.
  1178 +
  1179 + Examples
  1180 + ========
  1181 +
  1182 + >>> from sympy.stats import GammaInverse, density, cdf, E, variance
  1183 + >>> from sympy import Symbol, pprint
  1184 +
  1185 + >>> a = Symbol("a", positive=True)
  1186 + >>> b = Symbol("b", positive=True)
  1187 +
  1188 + >>> X = GammaInverse("x", a, b)
  1189 +
  1190 + >>> D = density(X)
  1191 + >>> pprint(D, use_unicode=False)
  1192 + / -b\
  1193 + | --|
  1194 + | -a - 1 a x |
  1195 + | x *b *e |
  1196 + Lambda|x, --------------|
  1197 + \ gamma(a) /
  1198 +
  1199 + References
  1200 + ==========
  1201 +
  1202 + [1] http://en.wikipedia.org/wiki/Inverse-gamma_distribution
  1203 + """
  1204 +
  1205 + return GammaInversePSpace(name, a, b).value
  1206 +
  1207 +#-------------------------------------------------------------------------------
  1208 +# Kumaraswamy distribution -----------------------------------------------------
  1209 +
  1210 +class KumaraswamyPSpace(SingleContinuousPSpace):
  1211 + def __new__(cls, name, a, b):
  1212 + a, b = sympify(a), sympify(b)
  1213 +
  1214 + _value_check(a > 0, "a must be positive")
  1215 + _value_check(b > 0, "b must be positive")
  1216 +
  1217 + x = Symbol(name)
  1218 + pdf = a * b * x**(a-1) * (1-x**a)**(b-1)
  1219 +
  1220 + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, 1))
  1221 + obj.a = a
  1222 + obj.b = b
  1223 + return obj
  1224 +
  1225 +def Kumaraswamy(name, a, b):
  1226 + r"""
  1227 + Create a Continuous Random Variable with a Kumaraswamy distribution.
  1228 +
  1229 + The density of the Kumaraswamy distribution is given by
  1230 +
  1231 + .. math::
  1232 + f(x) := a b x^{a-1} (1-x^a)^{b-1}
  1233 +
  1234 + with :math:`x \in [0,1]`.
  1235 +
  1236 + Parameters
  1237 + ==========
  1238 +
  1239 + a : Real number, `a` > 0 a shape
  1240 + b : Real number, `b` > 0 a shape
  1241 +
  1242 + Returns
  1243 + =======
  1244 +
  1245 + A RandomSymbol.
  1246 +
  1247 + Examples
  1248 + ========
  1249 +
  1250 + >>> from sympy.stats import Kumaraswamy, density, E, variance
  1251 + >>> from sympy import Symbol, simplify, pprint
  1252 +
  1253 + >>> a = Symbol("a", positive=True)
  1254 + >>> b = Symbol("b", positive=True)
  1255 +
  1256 + >>> X = Kumaraswamy("x", a, b)
  1257 +
  1258 + >>> D = density(X)
  1259 + >>> pprint(D, use_unicode=False)
  1260 + / b - 1\
  1261 + | a - 1 / a \ |
  1262 + Lambda\x, x *a*b*\- x + 1/ /
  1263 +
  1264 + >>> simplify(E(X, meijerg=True))
  1265 + gamma(1 + 1/a)*gamma(b + 1)/gamma(b + 1 + 1/a)
  1266 +
  1267 + References
  1268 + ==========
  1269 +
  1270 + [1] http://en.wikipedia.org/wiki/Kumaraswamy_distribution
  1271 + """
  1272 +
  1273 + return KumaraswamyPSpace(name, a, b).value
  1274 +
  1275 +#-------------------------------------------------------------------------------
734 1276 # Laplace distribution ---------------------------------------------------------
735 1277
736 1278
@@ -1221,6 +1763,148 @@ def Pareto(name, xm, alpha):
1221 1763 return ParetoPSpace(name, xm, alpha).value
1222 1764
1223 1765 #-------------------------------------------------------------------------------
  1766 +# QuadraticU distribution ------------------------------------------------------
  1767 +
  1768 +class QuadraticUPSpace(SingleContinuousPSpace):
  1769 + def __new__(cls, name, a, b):
  1770 + a, b = sympify(a), sympify(b)
  1771 + alpha = 12 / (b-a)**3
  1772 + beta = (a+b) / 2
  1773 +
  1774 + x = Symbol(name)
  1775 + pdf = Piecewise(
  1776 + (alpha * (x-beta)**2, And(a<=x, x<=b)),
  1777 + (S.Zero, True))
  1778 +
  1779 + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(a, b))
  1780 + obj.a = a
  1781 + obj.b = b
  1782 + return obj
  1783 +
  1784 +def QuadraticU(name, a, b):
  1785 + r"""
  1786 + Create a Continuous Random Variable with a U-quadratic distribution.
  1787 +
  1788 + The density of the U-quadratic distribution is given by
  1789 +
  1790 + .. math::
  1791 + f(x) := \alpha (x-\beta)^2
  1792 +
  1793 + with :math:`x \in [a,b]`.
  1794 +
  1795 + Parameters
  1796 + ==========
  1797 +
  1798 + a : Real number
  1799 + b : Real number, :math:`a < b`
  1800 +
  1801 + Returns
  1802 + =======
  1803 +
  1804 + A RandomSymbol.
  1805 +
  1806 + Examples
  1807 + ========
  1808 +
  1809 + >>> from sympy.stats import QuadraticU, density, E, variance
  1810 + >>> from sympy import Symbol, simplify, factor, pprint
  1811 +
  1812 + >>> a = Symbol("a", real=True)
  1813 + >>> b = Symbol("b", real=True)
  1814 +
  1815 + >>> X = QuadraticU("x", a, b)
  1816 +
  1817 + >>> D = density(X)
  1818 + >>> pprint(D, use_unicode=False)
  1819 + / / 2 \
  1820 + | | / a b\ |
  1821 + | |12*|x - - - -| |
  1822 + | | \ 2 2/ |
  1823 + Lambda|x, <--------------- for And(x <= b, a <= x)|
  1824 + | | 3 |
  1825 + | | (-a + b) |
  1826 + | | |
  1827 + \ \ 0 otherwise /
  1828 +
  1829 + References
  1830 + ==========
  1831 +
  1832 + [1] http://en.wikipedia.org/wiki/U-quadratic_distribution
  1833 + """
  1834 +
  1835 + return QuadraticUPSpace(name, a, b).value
  1836 +
  1837 +#-------------------------------------------------------------------------------
  1838 +# RaisedCosine distribution ----------------------------------------------------
  1839 +
  1840 +class RaisedCosinePSpace(SingleContinuousPSpace):
  1841 + def __new__(cls, name, mu, s):
  1842 + mu, s = sympify(mu), sympify(s)
  1843 +
  1844 + _value_check(s > 0, "s must be positive")
  1845 +
  1846 + x = Symbol(name)
  1847 + pdf = Piecewise(
  1848 + ((1+cos(pi*(x-mu)/s)) / (2*s), And(mu-s<=x, x<=mu+s)),
  1849 + (S.Zero, True))
  1850 +
  1851 + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(mu-s, mu+s))
  1852 + obj.mu = mu
  1853 + obj.s = s
  1854 + return obj
  1855 +
  1856 +def RaisedCosine(name, mu, s):
  1857 + r"""
  1858 + Create a Continuous Random Variable with a raised cosine distribution.
  1859 +
  1860 + The density of the raised cosine distribution is given by
  1861 +
  1862 + .. math::
  1863 + f(x) := \frac{1}{2s}\left(1+\cos\left(\frac{x-\mu}{s}\pi\right)\right)
  1864 +
  1865 + with :math:`x \in [\mu-s,\mu+s]`.
  1866 +
  1867 + Parameters
  1868 + ==========
  1869 +
  1870 + mu : Real number
  1871 + s : Real number, `s` > 0
  1872 +
  1873 + Returns
  1874 + =======
  1875 +
  1876 + A RandomSymbol.
  1877 +
  1878 + Examples
  1879 + ========
  1880 +
  1881 + >>> from sympy.stats import RaisedCosine, density, E, variance
  1882 + >>> from sympy import Symbol, simplify, pprint
  1883 +
  1884 + >>> mu = Symbol("mu", real=True)
  1885 + >>> s = Symbol("s", positive=True)
  1886 +
  1887 + >>> X = RaisedCosine("x", mu, s)
  1888 +
  1889 + >>> D = density(X)
  1890 + >>> pprint(D, use_unicode=False)
  1891 + / / /pi*(x - mu)\ \
  1892 + | |cos|-----------| + 1 |
  1893 + | | \ s / |
  1894 + Lambda|x, <-------------------- for And(x <= mu + s, mu - s <= x)|
  1895 + | | 2*s |
  1896 + | | |
  1897 + \ \ 0 otherwise /
  1898 +
  1899 + References
  1900 + ==========
  1901 +
  1902 + [1] http://en.wikipedia.org/wiki/Raised_cosine_distribution
  1903 + """
  1904 +
  1905 + return RaisedCosinePSpace(name, mu, s).value
  1906 +
  1907 +#-------------------------------------------------------------------------------
1224 1908 # Rayleigh distribution --------------------------------------------------------
1225 1909
1226 1910
@@ -1594,6 +2278,72 @@ def UniformSum(name, n):
1594 2278 return UniformSumPSpace(name, n).value
1595 2279
1596 2280 #-------------------------------------------------------------------------------
  2281 +# VonMises distribution --------------------------------------------------------
  2282 +
  2283 +class VonMisesPSpace(SingleContinuousPSpace):
  2284 + def __new__(cls, name, mu, k):
  2285 + mu, k = sympify(mu), sympify(k)
  2286 +
  2287 + _value_check(k > 0, "k must be positive")
  2288 +
  2289 + x = Symbol(name)
  2290 + pdf = exp(k*cos(x-mu)) / (2*pi*besseli(0, k))
  2291 +
  2292 + obj = SingleContinuousPSpace.__new__(cls, x, pdf, set=Interval(0, 2*pi))
  2293 + obj.mu = mu
  2294 + obj.k = k
  2295 + return obj
  2296 +
  2297 +def VonMises(name, mu, k):
  2298 + r"""
  2299 + Create a Continuous Random Variable with a von Mises distribution.
  2300 +
  2301 + The density of the von Mises distribution is given by
  2302 +
  2303 + .. math::
  2304 + f(x) := \frac{e^{\kappa\cos(x-\mu)}}{2\pi I_0(\kappa)}
  2305 +
  2306 + with :math:`x \in [0,2\pi]`.
  2307 +
  2308 + Parameters
  2309 + ==========
  2310 +
  2311 + mu : Real number, measure of location
  2312 + k : Real number, measure of concentration
  2313 +
  2314 + Returns
  2315 + =======
  2316 +
  2317 + A RandomSymbol.
  2318 +
  2319 + Examples
  2320 + ========
  2321 +
  2322 + >>> from sympy.stats import VonMises, density, E, variance
  2323 + >>> from sympy import Symbol, simplify, pprint
  2324 +
  2325 + >>> mu = Symbol("mu")
  2326 + >>> k = Symbol("k", positive=True)
  2327 +
  2328 + >