1from sympy import E as e
2from sympy import (Symbol, Abs, exp, expint, S, pi, simplify, Interval, erf, erfc, Ne,
3                   EulerGamma, Eq, log, lowergamma, uppergamma, symbols, sqrt, And,
4                   gamma, beta, Piecewise, Integral, sin, cos, tan, atan, sinh, cosh,
5                   besseli, floor, expand_func, Rational, I, re, Lambda, asin,
6                   im, lambdify, hyper, diff, Or, Mul, sign, Dummy, Sum,
7                   factorial, binomial, erfi, besselj, besselk)
8from sympy.functions.special.error_functions import erfinv
9from sympy.functions.special.hyper import meijerg
10from sympy.sets.sets import FiniteSet, Complement, Intersection
11from sympy.stats import (P, E, where, density, variance, covariance, skewness, kurtosis, median,
12                         given, pspace, cdf, characteristic_function, moment_generating_function,
13                         ContinuousRV, Arcsin, Benini, Beta, BetaNoncentral, BetaPrime,
14                         Cauchy, Chi, ChiSquared, ChiNoncentral, Dagum, Erlang, ExGaussian,
15                         Exponential, ExponentialPower, FDistribution, FisherZ, Frechet, Gamma,
16                         GammaInverse, Gompertz, Gumbel, Kumaraswamy, Laplace, Levy, Logistic, LogCauchy,
17                         LogLogistic, LogitNormal, LogNormal, Maxwell, Moyal, Nakagami, Normal, GaussianInverse,
18                         Pareto, PowerFunction, QuadraticU, RaisedCosine, Rayleigh, Reciprocal, ShiftedGompertz, StudentT,
19                         Trapezoidal, Triangular, Uniform, UniformSum, VonMises, Weibull, coskewness,
20                         WignerSemicircle, Wald, correlation, moment, cmoment, smoment, quantile,
21                         Lomax, BoundedPareto)
22
23from sympy.stats.crv_types import NormalDistribution, ExponentialDistribution, ContinuousDistributionHandmade
24from sympy.stats.joint_rv_types import MultivariateLaplaceDistribution, MultivariateNormalDistribution
25from sympy.stats.crv import SingleContinuousPSpace, SingleContinuousDomain
26from sympy.stats.compound_rv import CompoundPSpace
27from sympy.stats.symbolic_probability import Probability
28from sympy.testing.pytest import raises, XFAIL, slow, ignore_warnings
29from sympy.testing.randtest import verify_numerically as tn
30
31oo = S.Infinity
32
33x, y, z = map(Symbol, 'xyz')
34
35def test_single_normal():
36    mu = Symbol('mu', real=True)
37    sigma = Symbol('sigma', positive=True)
38    X = Normal('x', 0, 1)
39    Y = X*sigma + mu
40
41    assert E(Y) == mu
42    assert variance(Y) == sigma**2
43    pdf = density(Y)
44    x = Symbol('x', real=True)
45    assert (pdf(x) ==
46            2**S.Half*exp(-(x - mu)**2/(2*sigma**2))/(2*pi**S.Half*sigma))
47
48    assert P(X**2 < 1) == erf(2**S.Half/2)
49    ans = quantile(Y)(x)
50    assert ans == Complement(Intersection(FiniteSet(
51        sqrt(2)*sigma*(sqrt(2)*mu/(2*sigma)+ erfinv(2*x - 1))),
52        Interval(-oo, oo)), FiniteSet(mu))
53    assert E(X, Eq(X, mu)) == mu
54
55    assert median(X) == FiniteSet(0)
56    # issue 8248
57    assert X.pspace.compute_expectation(1).doit() == 1
58
59
60def test_conditional_1d():
61    X = Normal('x', 0, 1)
62    Y = given(X, X >= 0)
63    z = Symbol('z')
64
65    assert density(Y)(z) == 2 * density(X)(z)
66
67    assert Y.pspace.domain.set == Interval(0, oo)
68    assert E(Y) == sqrt(2) / sqrt(pi)
69
70    assert E(X**2) == E(Y**2)
71
72
73def test_ContinuousDomain():
74    X = Normal('x', 0, 1)
75    assert where(X**2 <= 1).set == Interval(-1, 1)
76    assert where(X**2 <= 1).symbol == X.symbol
77    where(And(X**2 <= 1, X >= 0)).set == Interval(0, 1)
78    raises(ValueError, lambda: where(sin(X) > 1))
79
80    Y = given(X, X >= 0)
81
82    assert Y.pspace.domain.set == Interval(0, oo)
83
84
85@slow
86def test_multiple_normal():
87    X, Y = Normal('x', 0, 1), Normal('y', 0, 1)
88    p = Symbol("p", positive=True)
89
90    assert E(X + Y) == 0
91    assert variance(X + Y) == 2
92    assert variance(X + X) == 4
93    assert covariance(X, Y) == 0
94    assert covariance(2*X + Y, -X) == -2*variance(X)
95    assert skewness(X) == 0
96    assert skewness(X + Y) == 0
97    assert kurtosis(X) == 3
98    assert kurtosis(X+Y) == 3
99    assert correlation(X, Y) == 0
100    assert correlation(X, X + Y) == correlation(X, X - Y)
101    assert moment(X, 2) == 1
102    assert cmoment(X, 3) == 0
103    assert moment(X + Y, 4) == 12
104    assert cmoment(X, 2) == variance(X)
105    assert smoment(X*X, 2) == 1
106    assert smoment(X + Y, 3) == skewness(X + Y)
107    assert smoment(X + Y, 4) == kurtosis(X + Y)
108    assert E(X, Eq(X + Y, 0)) == 0
109    assert variance(X, Eq(X + Y, 0)) == S.Half
110    assert quantile(X)(p) == sqrt(2)*erfinv(2*p - S.One)
111
112def test_symbolic():
113    mu1, mu2 = symbols('mu1 mu2', real=True)
114    s1, s2 = symbols('sigma1 sigma2', positive=True)
115    rate = Symbol('lambda', positive=True)
116    X = Normal('x', mu1, s1)
117    Y = Normal('y', mu2, s2)
118    Z = Exponential('z', rate)
119    a, b, c = symbols('a b c', real=True)
120
121    assert E(X) == mu1
122    assert E(X + Y) == mu1 + mu2
123    assert E(a*X + b) == a*E(X) + b
124    assert variance(X) == s1**2
125    assert variance(X + a*Y + b) == variance(X) + a**2*variance(Y)
126
127    assert E(Z) == 1/rate
128    assert E(a*Z + b) == a*E(Z) + b
129    assert E(X + a*Z + b) == mu1 + a/rate + b
130    assert median(X) == FiniteSet(mu1)
131
132
133def test_cdf():
134    X = Normal('x', 0, 1)
135
136    d = cdf(X)
137    assert P(X < 1) == d(1).rewrite(erfc)
138    assert d(0) == S.Half
139
140    d = cdf(X, X > 0)  # given X>0
141    assert d(0) == 0
142
143    Y = Exponential('y', 10)
144    d = cdf(Y)
145    assert d(-5) == 0
146    assert P(Y > 3) == 1 - d(3)
147
148    raises(ValueError, lambda: cdf(X + Y))
149
150    Z = Exponential('z', 1)
151    f = cdf(Z)
152    assert f(z) == Piecewise((1 - exp(-z), z >= 0), (0, True))
153
154
155def test_characteristic_function():
156    X = Uniform('x', 0, 1)
157
158    cf = characteristic_function(X)
159    assert cf(1) == -I*(-1 + exp(I))
160
161    Y = Normal('y', 1, 1)
162    cf = characteristic_function(Y)
163    assert cf(0) == 1
164    assert cf(1) == exp(I - S.Half)
165
166    Z = Exponential('z', 5)
167    cf = characteristic_function(Z)
168    assert cf(0) == 1
169    assert cf(1).expand() == Rational(25, 26) + I*Rational(5, 26)
170
171    X = GaussianInverse('x', 1, 1)
172    cf = characteristic_function(X)
173    assert cf(0) == 1
174    assert cf(1) == exp(1 - sqrt(1 - 2*I))
175
176    X = ExGaussian('x', 0, 1, 1)
177    cf = characteristic_function(X)
178    assert cf(0) == 1
179    assert cf(1) == (1 + I)*exp(Rational(-1, 2))/2
180
181    L = Levy('x', 0, 1)
182    cf = characteristic_function(L)
183    assert cf(0) == 1
184    assert cf(1) == exp(-sqrt(2)*sqrt(-I))
185
186
187def test_moment_generating_function():
188    t = symbols('t', positive=True)
189
190    # Symbolic tests
191    a, b, c = symbols('a b c')
192
193    mgf = moment_generating_function(Beta('x', a, b))(t)
194    assert mgf == hyper((a,), (a + b,), t)
195
196    mgf = moment_generating_function(Chi('x', a))(t)
197    assert mgf == sqrt(2)*t*gamma(a/2 + S.Half)*\
198        hyper((a/2 + S.Half,), (Rational(3, 2),), t**2/2)/gamma(a/2) +\
199        hyper((a/2,), (S.Half,), t**2/2)
200
201    mgf = moment_generating_function(ChiSquared('x', a))(t)
202    assert mgf == (1 - 2*t)**(-a/2)
203
204    mgf = moment_generating_function(Erlang('x', a, b))(t)
205    assert mgf == (1 - t/b)**(-a)
206
207    mgf = moment_generating_function(ExGaussian("x", a, b, c))(t)
208    assert mgf == exp(a*t + b**2*t**2/2)/(1 - t/c)
209
210    mgf = moment_generating_function(Exponential('x', a))(t)
211    assert mgf == a/(a - t)
212
213    mgf = moment_generating_function(Gamma('x', a, b))(t)
214    assert mgf == (-b*t + 1)**(-a)
215
216    mgf = moment_generating_function(Gumbel('x', a, b))(t)
217    assert mgf == exp(b*t)*gamma(-a*t + 1)
218
219    mgf = moment_generating_function(Gompertz('x', a, b))(t)
220    assert mgf == b*exp(b)*expint(t/a, b)
221
222    mgf = moment_generating_function(Laplace('x', a, b))(t)
223    assert mgf == exp(a*t)/(-b**2*t**2 + 1)
224
225    mgf = moment_generating_function(Logistic('x', a, b))(t)
226    assert mgf == exp(a*t)*beta(-b*t + 1, b*t + 1)
227
228    mgf = moment_generating_function(Normal('x', a, b))(t)
229    assert mgf == exp(a*t + b**2*t**2/2)
230
231    mgf = moment_generating_function(Pareto('x', a, b))(t)
232    assert mgf == b*(-a*t)**b*uppergamma(-b, -a*t)
233
234    mgf = moment_generating_function(QuadraticU('x', a, b))(t)
235    assert str(mgf) == ("(3*(t*(-4*b + (a + b)**2) + 4)*exp(b*t) - "
236    "3*(t*(a**2 + 2*a*(b - 2) + b**2) + 4)*exp(a*t))/(t**2*(a - b)**3)")
237
238    mgf = moment_generating_function(RaisedCosine('x', a, b))(t)
239    assert mgf == pi**2*exp(a*t)*sinh(b*t)/(b*t*(b**2*t**2 + pi**2))
240
241    mgf = moment_generating_function(Rayleigh('x', a))(t)
242    assert mgf == sqrt(2)*sqrt(pi)*a*t*(erf(sqrt(2)*a*t/2) + 1)\
243        *exp(a**2*t**2/2)/2 + 1
244
245    mgf = moment_generating_function(Triangular('x', a, b, c))(t)
246    assert str(mgf) == ("(-2*(-a + b)*exp(c*t) + 2*(-a + c)*exp(b*t) + "
247    "2*(b - c)*exp(a*t))/(t**2*(-a + b)*(-a + c)*(b - c))")
248
249    mgf = moment_generating_function(Uniform('x', a, b))(t)
250    assert mgf == (-exp(a*t) + exp(b*t))/(t*(-a + b))
251
252    mgf = moment_generating_function(UniformSum('x', a))(t)
253    assert mgf == ((exp(t) - 1)/t)**a
254
255    mgf = moment_generating_function(WignerSemicircle('x', a))(t)
256    assert mgf == 2*besseli(1, a*t)/(a*t)
257
258    # Numeric tests
259
260    mgf = moment_generating_function(Beta('x', 1, 1))(t)
261    assert mgf.diff(t).subs(t, 1) == hyper((2,), (3,), 1)/2
262
263    mgf = moment_generating_function(Chi('x', 1))(t)
264    assert mgf.diff(t).subs(t, 1) == sqrt(2)*hyper((1,), (Rational(3, 2),), S.Half
265    )/sqrt(pi) + hyper((Rational(3, 2),), (Rational(3, 2),), S.Half) + 2*sqrt(2)*hyper((2,),
266    (Rational(5, 2),), S.Half)/(3*sqrt(pi))
267
268    mgf = moment_generating_function(ChiSquared('x', 1))(t)
269    assert mgf.diff(t).subs(t, 1) == I
270
271    mgf = moment_generating_function(Erlang('x', 1, 1))(t)
272    assert mgf.diff(t).subs(t, 0) == 1
273
274    mgf = moment_generating_function(ExGaussian("x", 0, 1, 1))(t)
275    assert mgf.diff(t).subs(t, 2) == -exp(2)
276
277    mgf = moment_generating_function(Exponential('x', 1))(t)
278    assert mgf.diff(t).subs(t, 0) == 1
279
280    mgf = moment_generating_function(Gamma('x', 1, 1))(t)
281    assert mgf.diff(t).subs(t, 0) == 1
282
283    mgf = moment_generating_function(Gumbel('x', 1, 1))(t)
284    assert mgf.diff(t).subs(t, 0) == EulerGamma + 1
285
286    mgf = moment_generating_function(Gompertz('x', 1, 1))(t)
287    assert mgf.diff(t).subs(t, 1) == -e*meijerg(((), (1, 1)),
288    ((0, 0, 0), ()), 1)
289
290    mgf = moment_generating_function(Laplace('x', 1, 1))(t)
291    assert mgf.diff(t).subs(t, 0) == 1
292
293    mgf = moment_generating_function(Logistic('x', 1, 1))(t)
294    assert mgf.diff(t).subs(t, 0) == beta(1, 1)
295
296    mgf = moment_generating_function(Normal('x', 0, 1))(t)
297    assert mgf.diff(t).subs(t, 1) == exp(S.Half)
298
299    mgf = moment_generating_function(Pareto('x', 1, 1))(t)
300    assert mgf.diff(t).subs(t, 0) == expint(1, 0)
301
302    mgf = moment_generating_function(QuadraticU('x', 1, 2))(t)
303    assert mgf.diff(t).subs(t, 1) == -12*e - 3*exp(2)
304
305    mgf = moment_generating_function(RaisedCosine('x', 1, 1))(t)
306    assert mgf.diff(t).subs(t, 1) == -2*e*pi**2*sinh(1)/\
307    (1 + pi**2)**2 + e*pi**2*cosh(1)/(1 + pi**2)
308
309    mgf = moment_generating_function(Rayleigh('x', 1))(t)
310    assert mgf.diff(t).subs(t, 0) == sqrt(2)*sqrt(pi)/2
311
312    mgf = moment_generating_function(Triangular('x', 1, 3, 2))(t)
313    assert mgf.diff(t).subs(t, 1) == -e + exp(3)
314
315    mgf = moment_generating_function(Uniform('x', 0, 1))(t)
316    assert mgf.diff(t).subs(t, 1) == 1
317
318    mgf = moment_generating_function(UniformSum('x', 1))(t)
319    assert mgf.diff(t).subs(t, 1) == 1
320
321    mgf = moment_generating_function(WignerSemicircle('x', 1))(t)
322    assert mgf.diff(t).subs(t, 1) == -2*besseli(1, 1) + besseli(2, 1) +\
323        besseli(0, 1)
324
325
326def test_ContinuousRV():
327    pdf = sqrt(2)*exp(-x**2/2)/(2*sqrt(pi))  # Normal distribution
328    # X and Y should be equivalent
329    X = ContinuousRV(x, pdf, check=True)
330    Y = Normal('y', 0, 1)
331
332    assert variance(X) == variance(Y)
333    assert P(X > 0) == P(Y > 0)
334    Z = ContinuousRV(z, exp(-z), set=Interval(0, oo))
335    assert Z.pspace.domain.set == Interval(0, oo)
336    assert E(Z) == 1
337    assert P(Z > 5) == exp(-5)
338    raises(ValueError, lambda: ContinuousRV(z, exp(-z), set=Interval(0, 10), check=True))
339
340    # the correct pdf for Gamma(k, theta) but the integral in `check`
341    # integrates to something equivalent to 1 and not to 1 exactly
342    _x, k, theta = symbols("x k theta", positive=True)
343    pdf = 1/(gamma(k)*theta**k)*_x**(k-1)*exp(-_x/theta)
344    X = ContinuousRV(_x, pdf, set=Interval(0, oo))
345    Y = Gamma('y', k, theta)
346    assert (E(X) - E(Y)).simplify() == 0
347    assert (variance(X) - variance(Y)).simplify() == 0
348
349
350def test_arcsin():
351
352    a = Symbol("a", real=True)
353    b = Symbol("b", real=True)
354
355    X = Arcsin('x', a, b)
356    assert density(X)(x) == 1/(pi*sqrt((-x + b)*(x - a)))
357    assert cdf(X)(x) == Piecewise((0, a > x),
358                            (2*asin(sqrt((-a + x)/(-a + b)))/pi, b >= x),
359                            (1, True))
360    assert pspace(X).domain.set == Interval(a, b)
361
362def test_benini():
363    alpha = Symbol("alpha", positive=True)
364    beta = Symbol("beta", positive=True)
365    sigma = Symbol("sigma", positive=True)
366    X = Benini('x', alpha, beta, sigma)
367
368    assert density(X)(x) == ((alpha/x + 2*beta*log(x/sigma)/x)
369                          *exp(-alpha*log(x/sigma) - beta*log(x/sigma)**2))
370
371    assert pspace(X).domain.set == Interval(sigma, oo)
372    raises(NotImplementedError, lambda: moment_generating_function(X))
373    alpha = Symbol("alpha", nonpositive=True)
374    raises(ValueError, lambda: Benini('x', alpha, beta, sigma))
375
376    beta = Symbol("beta", nonpositive=True)
377    raises(ValueError, lambda: Benini('x', alpha, beta, sigma))
378
379    alpha = Symbol("alpha", positive=True)
380    raises(ValueError, lambda: Benini('x', alpha, beta, sigma))
381
382    beta = Symbol("beta", positive=True)
383    sigma = Symbol("sigma", nonpositive=True)
384    raises(ValueError, lambda: Benini('x', alpha, beta, sigma))
385
386def test_beta():
387    a, b = symbols('alpha beta', positive=True)
388    B = Beta('x', a, b)
389
390    assert pspace(B).domain.set == Interval(0, 1)
391    assert characteristic_function(B)(x) == hyper((a,), (a + b,), I*x)
392    assert density(B)(x) == x**(a - 1)*(1 - x)**(b - 1)/beta(a, b)
393
394    assert simplify(E(B)) == a / (a + b)
395    assert simplify(variance(B)) == a*b / (a**3 + 3*a**2*b + a**2 + 3*a*b**2 + 2*a*b + b**3 + b**2)
396
397    # Full symbolic solution is too much, test with numeric version
398    a, b = 1, 2
399    B = Beta('x', a, b)
400    assert expand_func(E(B)) == a / S(a + b)
401    assert expand_func(variance(B)) == (a*b) / S((a + b)**2 * (a + b + 1))
402    assert median(B) == FiniteSet(1 - 1/sqrt(2))
403
404def test_beta_noncentral():
405    a, b = symbols('a b', positive=True)
406    c = Symbol('c', nonnegative=True)
407    _k = Dummy('k')
408
409    X = BetaNoncentral('x', a, b, c)
410
411    assert pspace(X).domain.set == Interval(0, 1)
412
413    dens = density(X)
414    z = Symbol('z')
415
416    res = Sum( z**(_k + a - 1)*(c/2)**_k*(1 - z)**(b - 1)*exp(-c/2)/
417               (beta(_k + a, b)*factorial(_k)), (_k, 0, oo))
418    assert dens(z).dummy_eq(res)
419
420    # BetaCentral should not raise if the assumptions
421    # on the symbols can not be determined
422    a, b, c = symbols('a b c')
423    assert BetaNoncentral('x', a, b, c)
424
425    a = Symbol('a', positive=False, real=True)
426    raises(ValueError, lambda: BetaNoncentral('x', a, b, c))
427
428    a = Symbol('a', positive=True)
429    b = Symbol('b', positive=False, real=True)
430    raises(ValueError, lambda: BetaNoncentral('x', a, b, c))
431
432    a = Symbol('a', positive=True)
433    b = Symbol('b', positive=True)
434    c = Symbol('c', nonnegative=False, real=True)
435    raises(ValueError, lambda: BetaNoncentral('x', a, b, c))
436
437def test_betaprime():
438    alpha = Symbol("alpha", positive=True)
439
440    betap = Symbol("beta", positive=True)
441
442    X = BetaPrime('x', alpha, betap)
443    assert density(X)(x) == x**(alpha - 1)*(x + 1)**(-alpha - betap)/beta(alpha, betap)
444
445    alpha = Symbol("alpha", nonpositive=True)
446    raises(ValueError, lambda: BetaPrime('x', alpha, betap))
447
448    alpha = Symbol("alpha", positive=True)
449    betap = Symbol("beta", nonpositive=True)
450    raises(ValueError, lambda: BetaPrime('x', alpha, betap))
451    X = BetaPrime('x', 1, 1)
452    assert median(X) == FiniteSet(1)
453
454
455def test_BoundedPareto():
456    L, H = symbols('L, H', negative=True)
457    raises(ValueError, lambda: BoundedPareto('X', 1, L, H))
458    L, H = symbols('L, H', real=False)
459    raises(ValueError, lambda: BoundedPareto('X', 1, L, H))
460    L, H = symbols('L, H', positive=True)
461    raises(ValueError, lambda: BoundedPareto('X', -1, L, H))
462
463    X = BoundedPareto('X', 2, L, H)
464    assert X.pspace.domain.set == Interval(L, H)
465    assert density(X)(x) == 2*L**2/(x**3*(1 - L**2/H**2))
466    assert cdf(X)(x) == Piecewise((-H**2*L**2/(x**2*(H**2 - L**2)) \
467                            + H**2/(H**2 - L**2), L <= x), (0, True))
468    assert E(X).simplify() == 2*H*L/(H + L)
469    X = BoundedPareto('X', 1, 2, 4)
470    assert E(X).simplify() == log(16)
471    assert median(X) == FiniteSet(Rational(8, 3))
472    assert variance(X).simplify() == 8 - 16*log(2)**2
473
474
475def test_cauchy():
476    x0 = Symbol("x0", real=True)
477    gamma = Symbol("gamma", positive=True)
478    p = Symbol("p", positive=True)
479
480    X = Cauchy('x', x0, gamma)
481    # Tests the characteristic function
482    assert characteristic_function(X)(x) == exp(-gamma*Abs(x) + I*x*x0)
483    raises(NotImplementedError, lambda: moment_generating_function(X))
484    assert density(X)(x) == 1/(pi*gamma*(1 + (x - x0)**2/gamma**2))
485    assert diff(cdf(X)(x), x) == density(X)(x)
486    assert quantile(X)(p) == gamma*tan(pi*(p - S.Half)) + x0
487
488    x1 = Symbol("x1", real=False)
489    raises(ValueError, lambda: Cauchy('x', x1, gamma))
490    gamma = Symbol("gamma", nonpositive=True)
491    raises(ValueError, lambda: Cauchy('x', x0, gamma))
492    assert median(X) == FiniteSet(x0)
493
494def test_chi():
495    from sympy import I
496    k = Symbol("k", integer=True)
497
498    X = Chi('x', k)
499    assert density(X)(x) == 2**(-k/2 + 1)*x**(k - 1)*exp(-x**2/2)/gamma(k/2)
500
501    # Tests the characteristic function
502    assert characteristic_function(X)(x) == sqrt(2)*I*x*gamma(k/2 + S(1)/2)*hyper((k/2 + S(1)/2,),
503                                            (S(3)/2,), -x**2/2)/gamma(k/2) + hyper((k/2,), (S(1)/2,), -x**2/2)
504
505    # Tests the moment generating function
506    assert moment_generating_function(X)(x) == sqrt(2)*x*gamma(k/2 + S(1)/2)*hyper((k/2 + S(1)/2,),
507                                                (S(3)/2,), x**2/2)/gamma(k/2) + hyper((k/2,), (S(1)/2,), x**2/2)
508
509    k = Symbol("k", integer=True, positive=False)
510    raises(ValueError, lambda: Chi('x', k))
511
512    k = Symbol("k", integer=False, positive=True)
513    raises(ValueError, lambda: Chi('x', k))
514
515def test_chi_noncentral():
516    k = Symbol("k", integer=True)
517    l = Symbol("l")
518
519    X = ChiNoncentral("x", k, l)
520    assert density(X)(x) == (x**k*l*(x*l)**(-k/2)*
521                          exp(-x**2/2 - l**2/2)*besseli(k/2 - 1, x*l))
522
523    k = Symbol("k", integer=True, positive=False)
524    raises(ValueError, lambda: ChiNoncentral('x', k, l))
525
526    k = Symbol("k", integer=True, positive=True)
527    l = Symbol("l", nonpositive=True)
528    raises(ValueError, lambda: ChiNoncentral('x', k, l))
529
530    k = Symbol("k", integer=False)
531    l = Symbol("l", positive=True)
532    raises(ValueError, lambda: ChiNoncentral('x', k, l))
533
534
535def test_chi_squared():
536    k = Symbol("k", integer=True)
537    X = ChiSquared('x', k)
538
539    # Tests the characteristic function
540    assert characteristic_function(X)(x) == ((-2*I*x + 1)**(-k/2))
541
542    assert density(X)(x) == 2**(-k/2)*x**(k/2 - 1)*exp(-x/2)/gamma(k/2)
543    assert cdf(X)(x) == Piecewise((lowergamma(k/2, x/2)/gamma(k/2), x >= 0), (0, True))
544    assert E(X) == k
545    assert variance(X) == 2*k
546
547    X = ChiSquared('x', 15)
548    assert cdf(X)(3) == -14873*sqrt(6)*exp(Rational(-3, 2))/(5005*sqrt(pi)) + erf(sqrt(6)/2)
549
550    k = Symbol("k", integer=True, positive=False)
551    raises(ValueError, lambda: ChiSquared('x', k))
552
553    k = Symbol("k", integer=False, positive=True)
554    raises(ValueError, lambda: ChiSquared('x', k))
555
556
557def test_dagum():
558    p = Symbol("p", positive=True)
559    b = Symbol("b", positive=True)
560    a = Symbol("a", positive=True)
561
562    X = Dagum('x', p, a, b)
563    assert density(X)(x) == a*p*(x/b)**(a*p)*((x/b)**a + 1)**(-p - 1)/x
564    assert cdf(X)(x) == Piecewise(((1 + (x/b)**(-a))**(-p), x >= 0),
565                                    (0, True))
566
567    p = Symbol("p", nonpositive=True)
568    raises(ValueError, lambda: Dagum('x', p, a, b))
569
570    p = Symbol("p", positive=True)
571    b = Symbol("b", nonpositive=True)
572    raises(ValueError, lambda: Dagum('x', p, a, b))
573
574    b = Symbol("b", positive=True)
575    a = Symbol("a", nonpositive=True)
576    raises(ValueError, lambda: Dagum('x', p, a, b))
577    X = Dagum('x', 1 , 1, 1)
578    assert median(X) == FiniteSet(1)
579
580def test_erlang():
581    k = Symbol("k", integer=True, positive=True)
582    l = Symbol("l", positive=True)
583
584    X = Erlang("x", k, l)
585    assert density(X)(x) == x**(k - 1)*l**k*exp(-x*l)/gamma(k)
586    assert cdf(X)(x) == Piecewise((lowergamma(k, l*x)/gamma(k), x > 0),
587                               (0, True))
588
589
590def test_exgaussian():
591    m, z = symbols("m, z")
592    s, l = symbols("s, l", positive=True)
593    X = ExGaussian("x", m, s, l)
594
595    assert density(X)(z) == l*exp(l*(l*s**2 + 2*m - 2*z)/2) *\
596        erfc(sqrt(2)*(l*s**2 + m - z)/(2*s))/2
597
598    # Note: actual_output simplifies to expected_output.
599    # Ideally cdf(X)(z) would return expected_output
600    # expected_output = (erf(sqrt(2)*(l*s**2 + m - z)/(2*s)) - 1)*exp(l*(l*s**2 + 2*m - 2*z)/2)/2 - erf(sqrt(2)*(m - z)/(2*s))/2 + S.Half
601    u = l*(z - m)
602    v = l*s
603    GaussianCDF1 = cdf(Normal('x', 0, v))(u)
604    GaussianCDF2 = cdf(Normal('x', v**2, v))(u)
605    actual_output = GaussianCDF1 - exp(-u + (v**2/2) + log(GaussianCDF2))
606    assert cdf(X)(z) == actual_output
607    # assert simplify(actual_output) == expected_output
608
609    assert variance(X).expand() == s**2 + l**(-2)
610
611    assert skewness(X).expand() == 2/(l**3*s**2*sqrt(s**2 + l**(-2)) + l *
612                                      sqrt(s**2 + l**(-2)))
613
614
615def test_exponential():
616    rate = Symbol('lambda', positive=True)
617    X = Exponential('x', rate)
618    p = Symbol("p", positive=True, real=True, finite=True)
619
620    assert E(X) == 1/rate
621    assert variance(X) == 1/rate**2
622    assert skewness(X) == 2
623    assert skewness(X) == smoment(X, 3)
624    assert kurtosis(X) == 9
625    assert kurtosis(X) == smoment(X, 4)
626    assert smoment(2*X, 4) == smoment(X, 4)
627    assert moment(X, 3) == 3*2*1/rate**3
628    assert P(X > 0) is S.One
629    assert P(X > 1) == exp(-rate)
630    assert P(X > 10) == exp(-10*rate)
631    assert quantile(X)(p) == -log(1-p)/rate
632
633    assert where(X <= 1).set == Interval(0, 1)
634    Y = Exponential('y', 1)
635    assert median(Y) == FiniteSet(log(2))
636    #Test issue 9970
637    z = Dummy('z')
638    assert P(X > z) == exp(-z*rate)
639    assert P(X < z) == 0
640    #Test issue 10076 (Distribution with interval(0,oo))
641    x = Symbol('x')
642    _z = Dummy('_z')
643    b = SingleContinuousPSpace(x, ExponentialDistribution(2))
644
645    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
646        expected1 = Integral(2*exp(-2*_z), (_z, 3, oo))
647        assert b.probability(x > 3, evaluate=False).rewrite(Integral).dummy_eq(expected1)
648
649        expected2 = Integral(2*exp(-2*_z), (_z, 0, 4))
650        assert b.probability(x < 4, evaluate=False).rewrite(Integral).dummy_eq(expected2)
651    Y = Exponential('y', 2*rate)
652    assert coskewness(X, X, X) == skewness(X)
653    assert coskewness(X, Y + rate*X, Y + 2*rate*X) == \
654                        4/(sqrt(1 + 1/(4*rate**2))*sqrt(4 + 1/(4*rate**2)))
655    assert coskewness(X + 2*Y, Y + X, Y + 2*X, X > 3) == \
656                        sqrt(170)*Rational(9, 85)
657
658def test_exponential_power():
659    mu = Symbol('mu')
660    z = Symbol('z')
661    alpha = Symbol('alpha', positive=True)
662    beta = Symbol('beta', positive=True)
663
664    X = ExponentialPower('x', mu, alpha, beta)
665
666    assert density(X)(z) == beta*exp(-(Abs(mu - z)/alpha)
667                                     ** beta)/(2*alpha*gamma(1/beta))
668    assert cdf(X)(z) == S.Half + lowergamma(1/beta,
669                            (Abs(mu - z)/alpha)**beta)*sign(-mu + z)/\
670                                (2*gamma(1/beta))
671
672
673def test_f_distribution():
674    d1 = Symbol("d1", positive=True)
675    d2 = Symbol("d2", positive=True)
676
677    X = FDistribution("x", d1, d2)
678
679    assert density(X)(x) == (d2**(d2/2)*sqrt((d1*x)**d1*(d1*x + d2)**(-d1 - d2))
680                             /(x*beta(d1/2, d2/2)))
681
682    raises(NotImplementedError, lambda: moment_generating_function(X))
683    d1 = Symbol("d1", nonpositive=True)
684    raises(ValueError, lambda: FDistribution('x', d1, d1))
685
686    d1 = Symbol("d1", positive=True, integer=False)
687    raises(ValueError, lambda: FDistribution('x', d1, d1))
688
689    d1 = Symbol("d1", positive=True)
690    d2 = Symbol("d2", nonpositive=True)
691    raises(ValueError, lambda: FDistribution('x', d1, d2))
692
693    d2 = Symbol("d2", positive=True, integer=False)
694    raises(ValueError, lambda: FDistribution('x', d1, d2))
695
696
697def test_fisher_z():
698    d1 = Symbol("d1", positive=True)
699    d2 = Symbol("d2", positive=True)
700
701    X = FisherZ("x", d1, d2)
702    assert density(X)(x) == (2*d1**(d1/2)*d2**(d2/2)*(d1*exp(2*x) + d2)
703                             **(-d1/2 - d2/2)*exp(d1*x)/beta(d1/2, d2/2))
704
705def test_frechet():
706    a = Symbol("a", positive=True)
707    s = Symbol("s", positive=True)
708    m = Symbol("m", real=True)
709
710    X = Frechet("x", a, s=s, m=m)
711    assert density(X)(x) == a*((x - m)/s)**(-a - 1)*exp(-((x - m)/s)**(-a))/s
712    assert cdf(X)(x) == Piecewise((exp(-((-m + x)/s)**(-a)), m <= x), (0, True))
713
714@slow
715def test_gamma():
716    k = Symbol("k", positive=True)
717    theta = Symbol("theta", positive=True)
718
719    X = Gamma('x', k, theta)
720
721    # Tests characteristic function
722    assert characteristic_function(X)(x) == ((-I*theta*x + 1)**(-k))
723
724    assert density(X)(x) == x**(k - 1)*theta**(-k)*exp(-x/theta)/gamma(k)
725    assert cdf(X, meijerg=True)(z) == Piecewise(
726            (-k*lowergamma(k, 0)/gamma(k + 1) +
727                k*lowergamma(k, z/theta)/gamma(k + 1), z >= 0),
728            (0, True))
729
730    # assert simplify(variance(X)) == k*theta**2  # handled numerically below
731    assert E(X) == moment(X, 1)
732
733    k, theta = symbols('k theta', positive=True)
734    X = Gamma('x', k, theta)
735    assert E(X) == k*theta
736    assert variance(X) == k*theta**2
737    assert skewness(X).expand() == 2/sqrt(k)
738    assert kurtosis(X).expand() == 3 + 6/k
739
740    Y = Gamma('y', 2*k, 3*theta)
741    assert coskewness(X, theta*X + Y, k*X + Y).simplify() == \
742        2*531441**(-k)*sqrt(k)*theta*(3*3**(12*k) - 2*531441**k) \
743        /(sqrt(k**2 + 18)*sqrt(theta**2 + 18))
744
745def test_gamma_inverse():
746    a = Symbol("a", positive=True)
747    b = Symbol("b", positive=True)
748    X = GammaInverse("x", a, b)
749    assert density(X)(x) == x**(-a - 1)*b**a*exp(-b/x)/gamma(a)
750    assert cdf(X)(x) == Piecewise((uppergamma(a, b/x)/gamma(a), x > 0), (0, True))
751    assert characteristic_function(X)(x) == 2 * (-I*b*x)**(a/2) \
752            * besselk(a, 2*sqrt(b)*sqrt(-I*x))/gamma(a)
753    raises(NotImplementedError, lambda: moment_generating_function(X))
754
755def test_gompertz():
756    b = Symbol("b", positive=True)
757    eta = Symbol("eta", positive=True)
758
759    X = Gompertz("x", b, eta)
760
761    assert density(X)(x) == b*eta*exp(eta)*exp(b*x)*exp(-eta*exp(b*x))
762    assert cdf(X)(x) == 1 - exp(eta)*exp(-eta*exp(b*x))
763    assert diff(cdf(X)(x), x) == density(X)(x)
764
765
766def test_gumbel():
767    beta = Symbol("beta", positive=True)
768    mu = Symbol("mu")
769    x = Symbol("x")
770    y = Symbol("y")
771    X = Gumbel("x", beta, mu)
772    Y = Gumbel("y", beta, mu, minimum=True)
773    assert density(X)(x).expand() == \
774    exp(mu/beta)*exp(-x/beta)*exp(-exp(mu/beta)*exp(-x/beta))/beta
775    assert density(Y)(y).expand() == \
776    exp(-mu/beta)*exp(y/beta)*exp(-exp(-mu/beta)*exp(y/beta))/beta
777    assert cdf(X)(x).expand() == \
778    exp(-exp(mu/beta)*exp(-x/beta))
779    assert characteristic_function(X)(x) == exp(I*mu*x)*gamma(-I*beta*x + 1)
780
781def test_kumaraswamy():
782    a = Symbol("a", positive=True)
783    b = Symbol("b", positive=True)
784
785    X = Kumaraswamy("x", a, b)
786    assert density(X)(x) == x**(a - 1)*a*b*(-x**a + 1)**(b - 1)
787    assert cdf(X)(x) == Piecewise((0, x < 0),
788                                (-(-x**a + 1)**b + 1, x <= 1),
789                                (1, True))
790
791
792def test_laplace():
793    mu = Symbol("mu")
794    b = Symbol("b", positive=True)
795
796    X = Laplace('x', mu, b)
797
798    #Tests characteristic_function
799    assert characteristic_function(X)(x) == (exp(I*mu*x)/(b**2*x**2 + 1))
800
801    assert density(X)(x) == exp(-Abs(x - mu)/b)/(2*b)
802    assert cdf(X)(x) == Piecewise((exp((-mu + x)/b)/2, mu > x),
803                            (-exp((mu - x)/b)/2 + 1, True))
804    X = Laplace('x', [1, 2], [[1, 0], [0, 1]])
805    assert isinstance(pspace(X).distribution, MultivariateLaplaceDistribution)
806
807def test_levy():
808    mu = Symbol("mu", real=True)
809    c = Symbol("c", positive=True)
810
811    X = Levy('x', mu, c)
812    assert X.pspace.domain.set == Interval(mu, oo)
813    assert density(X)(x) == sqrt(c/(2*pi))*exp(-c/(2*(x - mu)))/((x - mu)**(S.One + S.Half))
814    assert cdf(X)(x) == erfc(sqrt(c/(2*(x - mu))))
815
816    raises(NotImplementedError, lambda: moment_generating_function(X))
817    mu = Symbol("mu", real=False)
818    raises(ValueError, lambda: Levy('x',mu,c))
819
820    c = Symbol("c", nonpositive=True)
821    raises(ValueError, lambda: Levy('x',mu,c))
822
823    mu = Symbol("mu", real=True)
824    raises(ValueError, lambda: Levy('x',mu,c))
825
826def test_logcauchy():
827    mu = Symbol("mu" , positive=True)
828    sigma = Symbol("sigma" , positive=True)
829
830    X = LogCauchy("x", mu, sigma)
831
832    assert density(X)(x) == sigma/(x*pi*(sigma**2 + (-mu + log(x))**2))
833    assert cdf(X)(x) == atan((log(x) - mu)/sigma)/pi + S.Half
834
835
836def test_logistic():
837    mu = Symbol("mu", real=True)
838    s = Symbol("s", positive=True)
839    p = Symbol("p", positive=True)
840
841    X = Logistic('x', mu, s)
842
843    #Tests characteristics_function
844    assert characteristic_function(X)(x) == \
845           (Piecewise((pi*s*x*exp(I*mu*x)/sinh(pi*s*x), Ne(x, 0)), (1, True)))
846
847    assert density(X)(x) == exp((-x + mu)/s)/(s*(exp((-x + mu)/s) + 1)**2)
848    assert cdf(X)(x) == 1/(exp((mu - x)/s) + 1)
849    assert quantile(X)(p) == mu - s*log(-S.One + 1/p)
850
851def test_loglogistic():
852    a, b = symbols('a b')
853    assert LogLogistic('x', a, b)
854
855    a = Symbol('a', negative=True)
856    b = Symbol('b', positive=True)
857    raises(ValueError, lambda: LogLogistic('x', a, b))
858
859    a = Symbol('a', positive=True)
860    b = Symbol('b', negative=True)
861    raises(ValueError, lambda: LogLogistic('x', a, b))
862
863    a, b, z, p = symbols('a b z p', positive=True)
864    X = LogLogistic('x', a, b)
865    assert density(X)(z) == b*(z/a)**(b - 1)/(a*((z/a)**b + 1)**2)
866    assert cdf(X)(z) == 1/(1 + (z/a)**(-b))
867    assert quantile(X)(p) == a*(p/(1 - p))**(1/b)
868
869    # Expectation
870    assert E(X) == Piecewise((S.NaN, b <= 1), (pi*a/(b*sin(pi/b)), True))
871    b = symbols('b', prime=True) # b > 1
872    X = LogLogistic('x', a, b)
873    assert E(X) == pi*a/(b*sin(pi/b))
874    X = LogLogistic('x', 1, 2)
875    assert median(X) == FiniteSet(1)
876
877def test_logitnormal():
878    mu = Symbol('mu', real=True)
879    s = Symbol('s', positive=True)
880    X = LogitNormal('x', mu, s)
881    x = Symbol('x')
882
883    assert density(X)(x) == sqrt(2)*exp(-(-mu + log(x/(1 - x)))**2/(2*s**2))/(2*sqrt(pi)*s*x*(1 - x))
884    assert cdf(X)(x) == erf(sqrt(2)*(-mu + log(x/(1 - x)))/(2*s))/2 + S(1)/2
885
886def test_lognormal():
887    mean = Symbol('mu', real=True)
888    std = Symbol('sigma', positive=True)
889    X = LogNormal('x', mean, std)
890    # The sympy integrator can't do this too well
891    #assert E(X) == exp(mean+std**2/2)
892    #assert variance(X) == (exp(std**2)-1) * exp(2*mean + std**2)
893
894    # The sympy integrator can't do this too well
895    #assert E(X) ==
896    raises(NotImplementedError, lambda: moment_generating_function(X))
897    mu = Symbol("mu", real=True)
898    sigma = Symbol("sigma", positive=True)
899
900    X = LogNormal('x', mu, sigma)
901    assert density(X)(x) == (sqrt(2)*exp(-(-mu + log(x))**2
902                                    /(2*sigma**2))/(2*x*sqrt(pi)*sigma))
903    # Tests cdf
904    assert cdf(X)(x) == Piecewise(
905                        (erf(sqrt(2)*(-mu + log(x))/(2*sigma))/2
906                        + S(1)/2, x > 0), (0, True))
907
908    X = LogNormal('x', 0, 1)  # Mean 0, standard deviation 1
909    assert density(X)(x) == sqrt(2)*exp(-log(x)**2/2)/(2*x*sqrt(pi))
910
911
912def test_Lomax():
913    a, l = symbols('a, l', negative=True)
914    raises(ValueError, lambda: Lomax('X', a , l))
915    a, l = symbols('a, l', real=False)
916    raises(ValueError, lambda: Lomax('X', a , l))
917
918    a, l = symbols('a, l', positive=True)
919    X = Lomax('X', a, l)
920    assert X.pspace.domain.set == Interval(0, oo)
921    assert density(X)(x) == a*(1 + x/l)**(-a - 1)/l
922    assert cdf(X)(x) == Piecewise((1 - (1 + x/l)**(-a), x >= 0), (0, True))
923    a = 3
924    X = Lomax('X', a, l)
925    assert E(X) == l/2
926    assert median(X) == FiniteSet(l*(-1 + 2**Rational(1, 3)))
927    assert variance(X) == 3*l**2/4
928
929
930def test_maxwell():
931    a = Symbol("a", positive=True)
932
933    X = Maxwell('x', a)
934
935    assert density(X)(x) == (sqrt(2)*x**2*exp(-x**2/(2*a**2))/
936        (sqrt(pi)*a**3))
937    assert E(X) == 2*sqrt(2)*a/sqrt(pi)
938    assert variance(X) == -8*a**2/pi + 3*a**2
939    assert cdf(X)(x) == erf(sqrt(2)*x/(2*a)) - sqrt(2)*x*exp(-x**2/(2*a**2))/(sqrt(pi)*a)
940    assert diff(cdf(X)(x), x) == density(X)(x)
941
942def test_Moyal():
943    mu = Symbol('mu',real=False)
944    sigma = Symbol('sigma', real=True, positive=True)
945    raises(ValueError, lambda: Moyal('M',mu, sigma))
946
947    mu = Symbol('mu', real=True)
948    sigma = Symbol('sigma', real=True, negative=True)
949    raises(ValueError, lambda: Moyal('M',mu, sigma))
950
951    sigma = Symbol('sigma', real=True, positive=True)
952    M = Moyal('M', mu, sigma)
953    assert density(M)(z) == sqrt(2)*exp(-exp((mu - z)/sigma)/2
954                        - (-mu + z)/(2*sigma))/(2*sqrt(pi)*sigma)
955    assert cdf(M)(z).simplify() == 1 - erf(sqrt(2)*exp((mu - z)/(2*sigma))/2)
956    assert characteristic_function(M)(z) == 2**(-I*sigma*z)*exp(I*mu*z) \
957                        *gamma(-I*sigma*z + Rational(1, 2))/sqrt(pi)
958    assert E(M) == mu + EulerGamma*sigma + sigma*log(2)
959    assert moment_generating_function(M)(z) == 2**(-sigma*z)*exp(mu*z) \
960                        *gamma(-sigma*z + Rational(1, 2))/sqrt(pi)
961
962
963def test_nakagami():
964    mu = Symbol("mu", positive=True)
965    omega = Symbol("omega", positive=True)
966
967    X = Nakagami('x', mu, omega)
968    assert density(X)(x) == (2*x**(2*mu - 1)*mu**mu*omega**(-mu)
969                                *exp(-x**2*mu/omega)/gamma(mu))
970    assert simplify(E(X)) == (sqrt(mu)*sqrt(omega)
971                                            *gamma(mu + S.Half)/gamma(mu + 1))
972    assert simplify(variance(X)) == (
973    omega - omega*gamma(mu + S.Half)**2/(gamma(mu)*gamma(mu + 1)))
974    assert cdf(X)(x) == Piecewise(
975                                (lowergamma(mu, mu*x**2/omega)/gamma(mu), x > 0),
976                                (0, True))
977    X = Nakagami('x',1 ,1)
978    assert median(X) == FiniteSet(sqrt(log(2)))
979
980def test_gaussian_inverse():
981    # test for symbolic parameters
982    a, b = symbols('a b')
983    assert GaussianInverse('x', a, b)
984
985    # Inverse Gaussian distribution is also known as Wald distribution
986    # `GaussianInverse` can also be referred by the name `Wald`
987    a, b, z = symbols('a b z')
988    X = Wald('x', a, b)
989    assert density(X)(z) == sqrt(2)*sqrt(b/z**3)*exp(-b*(-a + z)**2/(2*a**2*z))/(2*sqrt(pi))
990
991    a, b = symbols('a b', positive=True)
992    z = Symbol('z', positive=True)
993
994    X = GaussianInverse('x', a, b)
995    assert density(X)(z) == sqrt(2)*sqrt(b)*sqrt(z**(-3))*exp(-b*(-a + z)**2/(2*a**2*z))/(2*sqrt(pi))
996    assert E(X) == a
997    assert variance(X).expand() == a**3/b
998    assert cdf(X)(z) == (S.Half - erf(sqrt(2)*sqrt(b)*(1 + z/a)/(2*sqrt(z)))/2)*exp(2*b/a) +\
999         erf(sqrt(2)*sqrt(b)*(-1 + z/a)/(2*sqrt(z)))/2 + S.Half
1000
1001    a = symbols('a', nonpositive=True)
1002    raises(ValueError, lambda: GaussianInverse('x', a, b))
1003
1004    a = symbols('a', positive=True)
1005    b = symbols('b', nonpositive=True)
1006    raises(ValueError, lambda: GaussianInverse('x', a, b))
1007
1008def test_pareto():
1009    xm, beta = symbols('xm beta', positive=True)
1010    alpha = beta + 5
1011    X = Pareto('x', xm, alpha)
1012
1013    dens = density(X)
1014
1015    #Tests cdf function
1016    assert cdf(X)(x) == \
1017           Piecewise((-x**(-beta - 5)*xm**(beta + 5) + 1, x >= xm), (0, True))
1018
1019    #Tests characteristic_function
1020    assert characteristic_function(X)(x) == \
1021           ((-I*x*xm)**(beta + 5)*(beta + 5)*uppergamma(-beta - 5, -I*x*xm))
1022
1023    assert dens(x) == x**(-(alpha + 1))*xm**(alpha)*(alpha)
1024
1025    assert simplify(E(X)) == alpha*xm/(alpha-1)
1026
1027    # computation of taylor series for MGF still too slow
1028    #assert simplify(variance(X)) == xm**2*alpha / ((alpha-1)**2*(alpha-2))
1029
1030
1031def test_pareto_numeric():
1032    xm, beta = 3, 2
1033    alpha = beta + 5
1034    X = Pareto('x', xm, alpha)
1035
1036    assert E(X) == alpha*xm/S(alpha - 1)
1037    assert variance(X) == xm**2*alpha / S((alpha - 1)**2*(alpha - 2))
1038    assert median(X) == FiniteSet(3*2**Rational(1, 7))
1039    # Skewness tests too slow. Try shortcutting function?
1040
1041
1042def test_PowerFunction():
1043    alpha = Symbol("alpha", nonpositive=True)
1044    a, b = symbols('a, b', real=True)
1045    raises (ValueError, lambda: PowerFunction('x', alpha, a, b))
1046
1047    a, b = symbols('a, b', real=False)
1048    raises (ValueError, lambda: PowerFunction('x', alpha, a, b))
1049
1050    alpha = Symbol("alpha", positive=True)
1051    a, b = symbols('a, b', real=True)
1052    raises (ValueError, lambda: PowerFunction('x', alpha, 5, 2))
1053
1054    X = PowerFunction('X', 2, a, b)
1055    assert density(X)(z) == (-2*a + 2*z)/(-a + b)**2
1056    assert cdf(X)(z) == Piecewise((a**2/(a**2 - 2*a*b + b**2) -
1057        2*a*z/(a**2 - 2*a*b + b**2) + z**2/(a**2 - 2*a*b + b**2), a <= z), (0, True))
1058
1059    X = PowerFunction('X', 2, 0, 1)
1060    assert density(X)(z) == 2*z
1061    assert cdf(X)(z) == Piecewise((z**2, z >= 0), (0,True))
1062    assert E(X) == Rational(2,3)
1063    assert P(X < 0) == 0
1064    assert P(X < 1) == 1
1065    assert median(X) == FiniteSet(1/sqrt(2))
1066
1067def test_raised_cosine():
1068    mu = Symbol("mu", real=True)
1069    s = Symbol("s", positive=True)
1070
1071    X = RaisedCosine("x", mu, s)
1072
1073    assert pspace(X).domain.set == Interval(mu - s, mu + s)
1074    #Tests characteristics_function
1075    assert characteristic_function(X)(x) == \
1076           Piecewise((exp(-I*pi*mu/s)/2, Eq(x, -pi/s)), (exp(I*pi*mu/s)/2, Eq(x, pi/s)), (pi**2*exp(I*mu*x)*sin(s*x)/(s*x*(-s**2*x**2 + pi**2)), True))
1077
1078    assert density(X)(x) == (Piecewise(((cos(pi*(x - mu)/s) + 1)/(2*s),
1079                          And(x <= mu + s, mu - s <= x)), (0, True)))
1080
1081
1082def test_rayleigh():
1083    sigma = Symbol("sigma", positive=True)
1084
1085    X = Rayleigh('x', sigma)
1086
1087    #Tests characteristic_function
1088    assert characteristic_function(X)(x) == (-sqrt(2)*sqrt(pi)*sigma*x*(erfi(sqrt(2)*sigma*x/2) - I)*exp(-sigma**2*x**2/2)/2 + 1)
1089
1090    assert density(X)(x) ==  x*exp(-x**2/(2*sigma**2))/sigma**2
1091    assert E(X) == sqrt(2)*sqrt(pi)*sigma/2
1092    assert variance(X) == -pi*sigma**2/2 + 2*sigma**2
1093    assert cdf(X)(x) == 1 - exp(-x**2/(2*sigma**2))
1094    assert diff(cdf(X)(x), x) == density(X)(x)
1095
1096def test_reciprocal():
1097    a = Symbol("a", real=True)
1098    b = Symbol("b", real=True)
1099
1100    X = Reciprocal('x', a, b)
1101    assert density(X)(x) == 1/(x*(-log(a) + log(b)))
1102    assert cdf(X)(x) == Piecewise((log(a)/(log(a) - log(b)) - log(x)/(log(a) - log(b)), a <= x), (0, True))
1103    X = Reciprocal('x', 5, 30)
1104
1105    assert E(X) == 25/(log(30) - log(5))
1106    assert P(X < 4) == S.Zero
1107    assert P(X < 20) == log(20) / (log(30) - log(5)) - log(5) / (log(30) - log(5))
1108    assert cdf(X)(10) == log(10) / (log(30) - log(5)) - log(5) / (log(30) - log(5))
1109
1110    a = symbols('a', nonpositive=True)
1111    raises(ValueError, lambda: Reciprocal('x', a, b))
1112
1113    a = symbols('a', positive=True)
1114    b = symbols('b', positive=True)
1115    raises(ValueError, lambda: Reciprocal('x', a + b, a))
1116
1117def test_shiftedgompertz():
1118    b = Symbol("b", positive=True)
1119    eta = Symbol("eta", positive=True)
1120    X = ShiftedGompertz("x", b, eta)
1121    assert density(X)(x) == b*(eta*(1 - exp(-b*x)) + 1)*exp(-b*x)*exp(-eta*exp(-b*x))
1122
1123
1124def test_studentt():
1125    nu = Symbol("nu", positive=True)
1126
1127    X = StudentT('x', nu)
1128    assert density(X)(x) == (1 + x**2/nu)**(-nu/2 - S.Half)/(sqrt(nu)*beta(S.Half, nu/2))
1129    assert cdf(X)(x) == S.Half + x*gamma(nu/2 + S.Half)*hyper((S.Half, nu/2 + S.Half),
1130                                (Rational(3, 2),), -x**2/nu)/(sqrt(pi)*sqrt(nu)*gamma(nu/2))
1131    raises(NotImplementedError, lambda: moment_generating_function(X))
1132
1133def test_trapezoidal():
1134    a = Symbol("a", real=True)
1135    b = Symbol("b", real=True)
1136    c = Symbol("c", real=True)
1137    d = Symbol("d", real=True)
1138
1139    X = Trapezoidal('x', a, b, c, d)
1140    assert density(X)(x) == Piecewise(((-2*a + 2*x)/((-a + b)*(-a - b + c + d)), (a <= x) & (x < b)),
1141                                      (2/(-a - b + c + d), (b <= x) & (x < c)),
1142                                      ((2*d - 2*x)/((-c + d)*(-a - b + c + d)), (c <= x) & (x <= d)),
1143                                      (0, True))
1144
1145    X = Trapezoidal('x', 0, 1, 2, 3)
1146    assert E(X) == Rational(3, 2)
1147    assert variance(X) == Rational(5, 12)
1148    assert P(X < 2) == Rational(3, 4)
1149    assert median(X) == FiniteSet(Rational(3, 2))
1150
1151def test_triangular():
1152    a = Symbol("a")
1153    b = Symbol("b")
1154    c = Symbol("c")
1155
1156    X = Triangular('x', a, b, c)
1157    assert pspace(X).domain.set == Interval(a, b)
1158    assert str(density(X)(x)) == ("Piecewise(((-2*a + 2*x)/((-a + b)*(-a + c)), (a <= x) & (c > x)), "
1159    "(2/(-a + b), Eq(c, x)), ((2*b - 2*x)/((-a + b)*(b - c)), (b >= x) & (c < x)), (0, True))")
1160
1161    #Tests moment_generating_function
1162    assert moment_generating_function(X)(x).expand() == \
1163    ((-2*(-a + b)*exp(c*x) + 2*(-a + c)*exp(b*x) + 2*(b - c)*exp(a*x))/(x**2*(-a + b)*(-a + c)*(b - c))).expand()
1164    assert str(characteristic_function(X)(x)) == \
1165    '(2*(-a + b)*exp(I*c*x) - 2*(-a + c)*exp(I*b*x) - 2*(b - c)*exp(I*a*x))/(x**2*(-a + b)*(-a + c)*(b - c))'
1166
1167def test_quadratic_u():
1168    a = Symbol("a", real=True)
1169    b = Symbol("b", real=True)
1170
1171    X = QuadraticU("x", a, b)
1172    Y = QuadraticU("x", 1, 2)
1173
1174    assert pspace(X).domain.set == Interval(a, b)
1175    # Tests _moment_generating_function
1176    assert moment_generating_function(Y)(1)  == -15*exp(2) + 27*exp(1)
1177    assert moment_generating_function(Y)(2) == -9*exp(4)/2 + 21*exp(2)/2
1178
1179    assert characteristic_function(Y)(1) == 3*I*(-1 + 4*I)*exp(I*exp(2*I))
1180    assert density(X)(x) == (Piecewise((12*(x - a/2 - b/2)**2/(-a + b)**3,
1181                          And(x <= b, a <= x)), (0, True)))
1182
1183
1184def test_uniform():
1185    l = Symbol('l', real=True)
1186    w = Symbol('w', positive=True)
1187    X = Uniform('x', l, l + w)
1188
1189    assert E(X) == l + w/2
1190    assert variance(X).expand() == w**2/12
1191
1192    # With numbers all is well
1193    X = Uniform('x', 3, 5)
1194    assert P(X < 3) == 0 and P(X > 5) == 0
1195    assert P(X < 4) == P(X > 4) == S.Half
1196    assert median(X) == FiniteSet(4)
1197
1198    z = Symbol('z')
1199    p = density(X)(z)
1200    assert p.subs(z, 3.7) == S.Half
1201    assert p.subs(z, -1) == 0
1202    assert p.subs(z, 6) == 0
1203
1204    c = cdf(X)
1205    assert c(2) == 0 and c(3) == 0
1206    assert c(Rational(7, 2)) == Rational(1, 4)
1207    assert c(5) == 1 and c(6) == 1
1208
1209@XFAIL
1210def test_uniform_P():
1211    """ This stopped working because SingleContinuousPSpace.compute_density no
1212    longer calls integrate on a DiracDelta but rather just solves directly.
1213    integrate used to call UniformDistribution.expectation which special-cased
1214    subsed out the Min and Max terms that Uniform produces
1215
1216    I decided to regress on this class for general cleanliness (and I suspect
1217    speed) of the algorithm.
1218    """
1219    l = Symbol('l', real=True)
1220    w = Symbol('w', positive=True)
1221    X = Uniform('x', l, l + w)
1222    assert P(X < l) == 0 and P(X > l + w) == 0
1223
1224
1225def test_uniformsum():
1226    n = Symbol("n", integer=True)
1227    _k = Dummy("k")
1228    x = Symbol("x")
1229
1230    X = UniformSum('x', n)
1231    res = Sum((-1)**_k*(-_k + x)**(n - 1)*binomial(n, _k), (_k, 0, floor(x)))/factorial(n - 1)
1232    assert density(X)(x).dummy_eq(res)
1233
1234    #Tests set functions
1235    assert X.pspace.domain.set == Interval(0, n)
1236
1237    #Tests the characteristic_function
1238    assert characteristic_function(X)(x) == (-I*(exp(I*x) - 1)/x)**n
1239
1240    #Tests the moment_generating_function
1241    assert moment_generating_function(X)(x) == ((exp(x) - 1)/x)**n
1242
1243
1244def test_von_mises():
1245    mu = Symbol("mu")
1246    k = Symbol("k", positive=True)
1247
1248    X = VonMises("x", mu, k)
1249    assert density(X)(x) == exp(k*cos(x - mu))/(2*pi*besseli(0, k))
1250
1251
1252def test_weibull():
1253    a, b = symbols('a b', positive=True)
1254    # FIXME: simplify(E(X)) seems to hang without extended_positive=True
1255    # On a Linux machine this had a rapid memory leak...
1256    # a, b = symbols('a b', positive=True)
1257    X = Weibull('x', a, b)
1258
1259    assert E(X).expand() == a * gamma(1 + 1/b)
1260    assert variance(X).expand() == (a**2 * gamma(1 + 2/b) - E(X)**2).expand()
1261    assert simplify(skewness(X)) == (2*gamma(1 + 1/b)**3 - 3*gamma(1 + 1/b)*gamma(1 + 2/b) + gamma(1 + 3/b))/(-gamma(1 + 1/b)**2 + gamma(1 + 2/b))**Rational(3, 2)
1262    assert simplify(kurtosis(X)) == (-3*gamma(1 + 1/b)**4 +\
1263        6*gamma(1 + 1/b)**2*gamma(1 + 2/b) - 4*gamma(1 + 1/b)*gamma(1 + 3/b) + gamma(1 + 4/b))/(gamma(1 + 1/b)**2 - gamma(1 + 2/b))**2
1264
1265def test_weibull_numeric():
1266    # Test for integers and rationals
1267    a = 1
1268    bvals = [S.Half, 1, Rational(3, 2), 5]
1269    for b in bvals:
1270        X = Weibull('x', a, b)
1271        assert simplify(E(X)) == expand_func(a * gamma(1 + 1/S(b)))
1272        assert simplify(variance(X)) == simplify(
1273            a**2 * gamma(1 + 2/S(b)) - E(X)**2)
1274        # Not testing Skew... it's slow with int/frac values > 3/2
1275
1276
1277def test_wignersemicircle():
1278    R = Symbol("R", positive=True)
1279
1280    X = WignerSemicircle('x', R)
1281    assert pspace(X).domain.set == Interval(-R, R)
1282    assert density(X)(x) == 2*sqrt(-x**2 + R**2)/(pi*R**2)
1283    assert E(X) == 0
1284
1285
1286    #Tests ChiNoncentralDistribution
1287    assert characteristic_function(X)(x) == \
1288           Piecewise((2*besselj(1, R*x)/(R*x), Ne(x, 0)), (1, True))
1289
1290
1291def test_input_value_assertions():
1292    a, b = symbols('a b')
1293    p, q = symbols('p q', positive=True)
1294    m, n = symbols('m n', positive=False, real=True)
1295
1296    raises(ValueError, lambda: Normal('x', 3, 0))
1297    raises(ValueError, lambda: Normal('x', m, n))
1298    Normal('X', a, p)  # No error raised
1299    raises(ValueError, lambda: Exponential('x', m))
1300    Exponential('Ex', p)  # No error raised
1301    for fn in [Pareto, Weibull, Beta, Gamma]:
1302        raises(ValueError, lambda: fn('x', m, p))
1303        raises(ValueError, lambda: fn('x', p, n))
1304        fn('x', p, q)  # No error raised
1305
1306
1307def test_unevaluated():
1308    X = Normal('x', 0, 1)
1309    k = Dummy('k')
1310    expr1 = Integral(sqrt(2)*k*exp(-k**2/2)/(2*sqrt(pi)), (k, -oo, oo))
1311    expr2 = Integral(sqrt(2)*exp(-k**2/2)/(2*sqrt(pi)), (k, 0, oo))
1312    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
1313        assert E(X, evaluate=False).rewrite(Integral).dummy_eq(expr1)
1314        assert E(X + 1, evaluate=False).rewrite(Integral).dummy_eq(expr1 + 1)
1315        assert P(X > 0, evaluate=False).rewrite(Integral).dummy_eq(expr2)
1316
1317    assert P(X > 0, X**2 < 1) == S.Half
1318
1319
1320def test_probability_unevaluated():
1321    T = Normal('T', 30, 3)
1322    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
1323        assert type(P(T > 33, evaluate=False)) == Probability
1324
1325
1326def test_density_unevaluated():
1327    X = Normal('X', 0, 1)
1328    Y = Normal('Y', 0, 2)
1329    assert isinstance(density(X+Y, evaluate=False)(z), Integral)
1330
1331
1332def test_NormalDistribution():
1333    nd = NormalDistribution(0, 1)
1334    x = Symbol('x')
1335    assert nd.cdf(x) == erf(sqrt(2)*x/2)/2 + S.Half
1336    assert nd.expectation(1, x) == 1
1337    assert nd.expectation(x, x) == 0
1338    assert nd.expectation(x**2, x) == 1
1339    #Test issue 10076
1340    a = SingleContinuousPSpace(x, NormalDistribution(2, 4))
1341    _z = Dummy('_z')
1342
1343    expected1 = Integral(sqrt(2)*exp(-(_z - 2)**2/32)/(8*sqrt(pi)),(_z, -oo, 1))
1344    assert a.probability(x < 1, evaluate=False).dummy_eq(expected1) is True
1345
1346    expected2 = Integral(sqrt(2)*exp(-(_z - 2)**2/32)/(8*sqrt(pi)),(_z, 1, oo))
1347    assert a.probability(x > 1, evaluate=False).dummy_eq(expected2) is True
1348
1349    b = SingleContinuousPSpace(x, NormalDistribution(1, 9))
1350
1351    expected3 = Integral(sqrt(2)*exp(-(_z - 1)**2/162)/(18*sqrt(pi)),(_z, 6, oo))
1352    assert b.probability(x > 6, evaluate=False).dummy_eq(expected3) is True
1353
1354    expected4 = Integral(sqrt(2)*exp(-(_z - 1)**2/162)/(18*sqrt(pi)),(_z, -oo, 6))
1355    assert b.probability(x < 6, evaluate=False).dummy_eq(expected4) is True
1356
1357
1358def test_random_parameters():
1359    mu = Normal('mu', 2, 3)
1360    meas = Normal('T', mu, 1)
1361    assert density(meas, evaluate=False)(z)
1362    assert isinstance(pspace(meas), CompoundPSpace)
1363    X = Normal('x', [1, 2], [[1, 0], [0, 1]])
1364    assert isinstance(pspace(X).distribution, MultivariateNormalDistribution)
1365    assert density(meas)(z).simplify() == sqrt(5)*exp(-z**2/20 + z/5 - S(1)/5)/(10*sqrt(pi))
1366
1367
1368def test_random_parameters_given():
1369    mu = Normal('mu', 2, 3)
1370    meas = Normal('T', mu, 1)
1371    assert given(meas, Eq(mu, 5)) == Normal('T', 5, 1)
1372
1373
1374def test_conjugate_priors():
1375    mu = Normal('mu', 2, 3)
1376    x = Normal('x', mu, 1)
1377    assert isinstance(simplify(density(mu, Eq(x, y), evaluate=False)(z)),
1378            Mul)
1379
1380
1381def test_difficult_univariate():
1382    """ Since using solve in place of deltaintegrate we're able to perform
1383    substantially more complex density computations on single continuous random
1384    variables """
1385    x = Normal('x', 0, 1)
1386    assert density(x**3)
1387    assert density(exp(x**2))
1388    assert density(log(x))
1389
1390
1391def test_issue_10003():
1392    X = Exponential('x', 3)
1393    G = Gamma('g', 1, 2)
1394    assert P(X < -1) is S.Zero
1395    assert P(G < -1) is S.Zero
1396
1397
1398@slow
1399def test_precomputed_cdf():
1400    x = symbols("x", real=True)
1401    mu = symbols("mu", real=True)
1402    sigma, xm, alpha = symbols("sigma xm alpha", positive=True)
1403    n = symbols("n", integer=True, positive=True)
1404    distribs = [
1405            Normal("X", mu, sigma),
1406            Pareto("P", xm, alpha),
1407            ChiSquared("C", n),
1408            Exponential("E", sigma),
1409            # LogNormal("L", mu, sigma),
1410    ]
1411    for X in distribs:
1412        compdiff = cdf(X)(x) - simplify(X.pspace.density.compute_cdf()(x))
1413        compdiff = simplify(compdiff.rewrite(erfc))
1414        assert compdiff == 0
1415
1416
1417@slow
1418def test_precomputed_characteristic_functions():
1419    import mpmath
1420
1421    def test_cf(dist, support_lower_limit, support_upper_limit):
1422        pdf = density(dist)
1423        t = Symbol('t')
1424
1425        # first function is the hardcoded CF of the distribution
1426        cf1 = lambdify([t], characteristic_function(dist)(t), 'mpmath')
1427
1428        # second function is the Fourier transform of the density function
1429        f = lambdify([x, t], pdf(x)*exp(I*x*t), 'mpmath')
1430        cf2 = lambda t: mpmath.quad(lambda x: f(x, t), [support_lower_limit, support_upper_limit], maxdegree=10)
1431
1432        # compare the two functions at various points
1433        for test_point in [2, 5, 8, 11]:
1434            n1 = cf1(test_point)
1435            n2 = cf2(test_point)
1436
1437            assert abs(re(n1) - re(n2)) < 1e-12
1438            assert abs(im(n1) - im(n2)) < 1e-12
1439
1440    test_cf(Beta('b', 1, 2), 0, 1)
1441    test_cf(Chi('c', 3), 0, mpmath.inf)
1442    test_cf(ChiSquared('c', 2), 0, mpmath.inf)
1443    test_cf(Exponential('e', 6), 0, mpmath.inf)
1444    test_cf(Logistic('l', 1, 2), -mpmath.inf, mpmath.inf)
1445    test_cf(Normal('n', -1, 5), -mpmath.inf, mpmath.inf)
1446    test_cf(RaisedCosine('r', 3, 1), 2, 4)
1447    test_cf(Rayleigh('r', 0.5), 0, mpmath.inf)
1448    test_cf(Uniform('u', -1, 1), -1, 1)
1449    test_cf(WignerSemicircle('w', 3), -3, 3)
1450
1451
1452def test_long_precomputed_cdf():
1453    x = symbols("x", real=True)
1454    distribs = [
1455            Arcsin("A", -5, 9),
1456            Dagum("D", 4, 10, 3),
1457            Erlang("E", 14, 5),
1458            Frechet("F", 2, 6, -3),
1459            Gamma("G", 2, 7),
1460            GammaInverse("GI", 3, 5),
1461            Kumaraswamy("K", 6, 8),
1462            Laplace("LA", -5, 4),
1463            Logistic("L", -6, 7),
1464            Nakagami("N", 2, 7),
1465            StudentT("S", 4)
1466            ]
1467    for distr in distribs:
1468        for _ in range(5):
1469            assert tn(diff(cdf(distr)(x), x), density(distr)(x), x, a=0, b=0, c=1, d=0)
1470
1471    US = UniformSum("US", 5)
1472    pdf01 = density(US)(x).subs(floor(x), 0).doit()   # pdf on (0, 1)
1473    cdf01 = cdf(US, evaluate=False)(x).subs(floor(x), 0).doit()   # cdf on (0, 1)
1474    assert tn(diff(cdf01, x), pdf01, x, a=0, b=0, c=1, d=0)
1475
1476
1477def test_issue_13324():
1478    X = Uniform('X', 0, 1)
1479    assert E(X, X > S.Half) == Rational(3, 4)
1480    assert E(X, X > 0) == S.Half
1481
1482def test_issue_20756():
1483    X = Uniform('X', -1, +1)
1484    Y = Uniform('Y', -1, +1)
1485    assert E(X * Y) == S.Zero
1486    assert E(X * ((Y + 1) - 1)) == S.Zero
1487    assert E(Y * (X*(X + 1) - X*X)) == S.Zero
1488
1489def test_FiniteSet_prob():
1490    E = Exponential('E', 3)
1491    N = Normal('N', 5, 7)
1492    assert P(Eq(E, 1)) is S.Zero
1493    assert P(Eq(N, 2)) is S.Zero
1494    assert P(Eq(N, x)) is S.Zero
1495
1496def test_prob_neq():
1497    E = Exponential('E', 4)
1498    X = ChiSquared('X', 4)
1499    assert P(Ne(E, 2)) == 1
1500    assert P(Ne(X, 4)) == 1
1501    assert P(Ne(X, 4)) == 1
1502    assert P(Ne(X, 5)) == 1
1503    assert P(Ne(E, x)) == 1
1504
1505def test_union():
1506    N = Normal('N', 3, 2)
1507    assert simplify(P(N**2 - N > 2)) == \
1508        -erf(sqrt(2))/2 - erfc(sqrt(2)/4)/2 + Rational(3, 2)
1509    assert simplify(P(N**2 - 4 > 0)) == \
1510        -erf(5*sqrt(2)/4)/2 - erfc(sqrt(2)/4)/2 + Rational(3, 2)
1511
1512def test_Or():
1513    N = Normal('N', 0, 1)
1514    assert simplify(P(Or(N > 2, N < 1))) == \
1515        -erf(sqrt(2))/2 - erfc(sqrt(2)/2)/2 + Rational(3, 2)
1516    assert P(Or(N < 0, N < 1)) == P(N < 1)
1517    assert P(Or(N > 0, N < 0)) == 1
1518
1519
1520def test_conditional_eq():
1521    E = Exponential('E', 1)
1522    assert P(Eq(E, 1), Eq(E, 1)) == 1
1523    assert P(Eq(E, 1), Eq(E, 2)) == 0
1524    assert P(E > 1, Eq(E, 2)) == 1
1525    assert P(E < 1, Eq(E, 2)) == 0
1526
1527def test_ContinuousDistributionHandmade():
1528    x = Symbol('x')
1529    z = Dummy('z')
1530    dens = Lambda(x, Piecewise((S.Half, (0<=x)&(x<1)), (0, (x>=1)&(x<2)),
1531        (S.Half, (x>=2)&(x<3)), (0, True)))
1532    dens = ContinuousDistributionHandmade(dens, set=Interval(0, 3))
1533    space = SingleContinuousPSpace(z, dens)
1534    assert dens.pdf == Lambda(x, Piecewise((1/2, (x >= 0) & (x < 1)),
1535        (0, (x >= 1) & (x < 2)), (1/2, (x >= 2) & (x < 3)), (0, True)))
1536    assert median(space.value) == Interval(1, 2)
1537    assert E(space.value) == Rational(3, 2)
1538    assert variance(space.value) == Rational(13, 12)
1539
1540
1541def test_issue_16318():
1542    # test compute_expectation function of the SingleContinuousDomain
1543    N = SingleContinuousDomain(x, Interval(0, 1))
1544    raises(ValueError, lambda: SingleContinuousDomain.compute_expectation(N, x+1, {x, y}))
1545