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