[docs]classContinuousDomain(RandomDomain):""" A domain with continuous support Represented using symbols and Intervals. """is_Continuous=Truedefas_boolean(self):raiseNotImplementedError("Not Implemented for generic Domains")

classSingleContinuousDomain(ContinuousDomain,SingleDomain):""" A univariate domain with continuous support Represented using a single symbol and interval. """def__new__(cls,symbol,set):assertsymbol.is_Symbolsymbols=FiniteSet(symbol)returnRandomDomain.__new__(cls,symbols,set)defintegrate(self,expr,variables=None,**kwargs):ifvariablesisNone:variables=self.symbolsifnotvariables:returnexprassertfrozenset(variables)==frozenset(self.symbols)# assumes only intervalsevaluate=kwargs.pop('evaluate',True)ifevaluate:returnintegrate(expr,(self.symbol,self.set),**kwargs)else:returnIntegral(expr,(self.symbol,self.set),**kwargs)defas_boolean(self):returnself.set.as_relational(self.symbol)classProductContinuousDomain(ProductDomain,ContinuousDomain):""" A collection of independent domains with continuous support """defintegrate(self,expr,variables=None,**kwargs):ifvariablesisNone:variables=self.symbolsfordomaininself.domains:domain_vars=frozenset(variables)&frozenset(domain.symbols)ifdomain_vars:expr=domain.integrate(expr,domain_vars,**kwargs)returnexprdefas_boolean(self):returnAnd(*[domain.as_boolean()fordomaininself.domains])classConditionalContinuousDomain(ContinuousDomain,ConditionalDomain):""" A domain with continuous support that has been further restricted by a condition such as x > 3 """defintegrate(self,expr,variables=None,**kwargs):ifvariablesisNone:variables=self.symbolsifnotvariables:returnexpr# Extract the full integralfullintgrl=self.fulldomain.integrate(expr,variables,evaluate=False)# separate into integrand and limitsintegrand,limits=fullintgrl.function,list(fullintgrl.limits)conditions=[self.condition]whileconditions:cond=conditions.pop()ifcond.is_Boolean:ifisinstance(cond,And):conditions.extend(cond.args)elifisinstance(cond,Or):raiseNotImplementedError("Or not implemented here")elifcond.is_Relational:ifcond.is_Equality:# Add the appropriate Delta to the integrandintegrand*=DiracDelta(cond.lhs-cond.rhs)else:symbols=cond.free_symbols&set(self.symbols)iflen(symbols)!=1:# Can't handle x > yraiseNotImplementedError("Multivariate Inequalities not yet implemented")# Can handle x > 0symbol=symbols.pop()# Find the limit with x, such as (x, -oo, oo)fori,limitinenumerate(limits):iflimit[0]==symbol:# Make condition into an Interval like [0, oo]cintvl=reduce_poly_inequalities_wrap(cond,symbol)# Make limit into an Interval like [-oo, oo]lintvl=Interval(limit[1],limit[2])# Intersect them to get [0, oo]intvl=cintvl.intersect(lintvl)# Put back into limits listlimits[i]=(symbol,intvl.left,intvl.right)else:raiseTypeError("Condition %s is not a relational or Boolean"%cond)evaluate=kwargs.pop('evaluate',True)ifevaluate:returnintegrate(integrand,*limits,**kwargs)returnIntegral(integrand,*limits,**kwargs)defas_boolean(self):returnAnd(self.fulldomain.as_boolean(),self.condition)@propertydefset(self):iflen(self.symbols)==1:return(self.fulldomain.set&reduce_poly_inequalities_wrap(self.condition,tuple(self.symbols)[0]))else:raiseNotImplementedError("Set of Conditional Domain not Implemented")

[docs]classContinuousPSpace(PSpace):""" A Continuous Probability Space Represents the likelihood of an event space defined over a continuum. Represented with a set of symbols and a probability density function. """is_Continuous=Truedefintegrate(self,expr,rvs=None,**kwargs):ifrvsisNone:rvs=self.valueselse:rvs=frozenset(rvs)expr=expr.subs(dict((rv,rv.symbol)forrvinrvs))domain_symbols=frozenset(rv.symbolforrvinrvs)returnself.domain.integrate(self.density*expr,domain_symbols,**kwargs)defcompute_density(self,expr,**kwargs):# Common case Density(X) where X in self.valuesifexprinself.values:# Marginalize all other random symbols out of the densitydensity=self.domain.integrate(self.density,set(rs.symbolforrsinself.values-frozenset((expr,))),**kwargs)returnLambda(expr.symbol,density)z=Dummy('z',real=True,bounded=True)returnLambda(z,self.integrate(DiracDelta(expr-z),**kwargs))@cacheitdefcompute_cdf(self,expr,**kwargs):ifnotself.domain.set.is_Interval:raiseValueError("CDF not well defined on multivariate expressions")d=self.compute_density(expr,**kwargs)x,z=symbols('x, z',real=True,bounded=True,cls=Dummy)left_bound=self.domain.set.start# CDF is integral of PDF from left bound to zcdf=integrate(d(x),(x,left_bound,z),**kwargs)# CDF Ensure that CDF left of left_bound is zerocdf=Piecewise((cdf,z>=left_bound),(0,True))returnLambda(z,cdf)defprobability(self,condition,**kwargs):z=Dummy('z',real=True,bounded=True)# Univariate case can be handled by wheretry:domain=self.where(condition)rv=[rvforrvinself.valuesifrv.symbol==domain.symbol][0]# Integrate out all other random variablespdf=self.compute_density(rv,**kwargs)# Integrate out the last variable over the special domainevaluate=kwargs.pop("evaluate",True)ifevaluate:returnintegrate(pdf(z),(z,domain.set),**kwargs)else:returnIntegral(pdf(z),(z,domain.set),**kwargs)# Other cases can be turned into univariate case# by computing a density handled by density computationexceptNotImplementedError:expr=condition.lhs-condition.rhsdensity=self.compute_density(expr,**kwargs)# Turn problem into univariate casespace=SingleContinuousPSpace(z,density(z))returnspace.probability(condition.__class__(space.value,0))defwhere(self,condition):rvs=frozenset(random_symbols(condition))ifnot(len(rvs)==1andrvs.issubset(self.values)):raiseNotImplementedError("Multiple continuous random variables not supported")rv=tuple(rvs)[0]interval=reduce_poly_inequalities_wrap(condition,rv)interval=interval.intersect(self.domain.set)returnSingleContinuousDomain(rv.symbol,interval)defconditional_space(self,condition,normalize=True,**kwargs):condition=condition.subs(dict((rv,rv.symbol)forrvinself.values))domain=ConditionalContinuousDomain(self.domain,condition)density=self.densityifnormalize:density=density/domain.integrate(density,**kwargs)returnContinuousPSpace(domain,density)

classSingleContinuousPSpace(ContinuousPSpace,SinglePSpace):""" A continuous probability space over a single univariate domain This class is commonly implemented by the various ContinuousRV types such as Normal, Exponential, Uniform, etc.... """def__new__(cls,symbol,density,set=Interval(-oo,oo)):domain=SingleContinuousDomain(sympify(symbol),set)obj=ContinuousPSpace.__new__(cls,domain,density)obj._cdf=Nonereturnobj@cacheitdef_inverse_cdf_expression(self):""" Inverse of the CDF See Also ======== compute_cdf sample """d=self.compute_cdf(self.value)x,z=symbols('x, z',real=True,positive=True,cls=Dummy)# Invert CDFtry:inverse_cdf=solve(d(x)-z,x)exceptNotImplementedError:inverse_cdf=Noneifnotinverse_cdforlen(inverse_cdf)!=1:raiseNotImplementedError("Could not invert CDF")returnLambda(z,inverse_cdf[0])defsample(self):""" Internal sample method Returns dictionary mapping RandomSymbol to realization value. """icdf=self._inverse_cdf_expression()return{self.value:icdf(random.uniform(0,1))}classProductContinuousPSpace(ProductPSpace,ContinuousPSpace):""" A collection of independent continuous probability spaces """@propertydefdensity(self):returnMul(*[space.densityforspaceinself.spaces])def_reduce_inequalities(conditions,var,**kwargs):try:returnreduce_poly_inequalities(conditions,var,**kwargs)exceptPolynomialError:raiseValueError("Reduction of condition failed %s\n"%conditions[0])defreduce_poly_inequalities_wrap(condition,var):ifcondition.is_Relational:return_reduce_inequalities([[condition]],var,relational=False)ifcondition.__class__isOr:return_reduce_inequalities([list(condition.args)],var,relational=False)ifcondition.__class__isAnd:intervals=[_reduce_inequalities([[arg]],var,relational=False)forargincondition.args]I=intervals[0]foriinintervals:I=I.intersect(i)returnI