From 1b7efdd676e2609abed1389e4b5df38a8358b8bd Mon Sep 17 00:00:00 2001
From: Stepan Roucka <stepan@roucka.eu>
Date: Sat, 30 Aug 2008 11:13:48 +0200
Subject: [PATCH] draft implementation of arbitrary PDF

Signed-off-by: Kirill Smelkov <kirr@landau.phys.spbu.ru>
---
 sympy/printing/pretty/pretty.py   |  17 +++++++
 sympy/printing/str.py             |   5 ++
 sympy/statistics/distributions.py | 100 ++++++++++++++++++++++++++++++++++++++
 3 files changed, 122 insertions(+)

diff --git a/sympy/printing/pretty/pretty.py b/sympy/printing/pretty/pretty.py
index 317f9a8..33e5bd8 100644
--- a/sympy/printing/pretty/pretty.py
+++ b/sympy/printing/pretty/pretty.py
@@ -110,6 +110,23 @@ class PrettyPrinter(Printer):
         pform = prettyForm(*stringPict.next(pform, f))
         return pform
 
+    def _print_PDF(self, pdf):
+        lim = self._print(pdf.pdf.args[0])
+        lim = prettyForm(*lim.right(', '))
+        lim = prettyForm(*lim.right(self._print(pdf.domain[0])))
+        lim = prettyForm(*lim.right(', '))
+        lim = prettyForm(*lim.right(self._print(pdf.domain[1])))
+        lim = prettyForm(*lim.parens())
+
+        f = self._print(pdf.pdf.args[1])
+        f = prettyForm(*f.right(', '))
+        f = prettyForm(*f.right(lim))
+        f = prettyForm(*f.parens())
+
+        pform = prettyForm('PDF')
+        pform = prettyForm(*pform.right(f))
+        return pform
+
     def _print_Integral(self, integral):
         f   = integral.function
 
diff --git a/sympy/printing/str.py b/sympy/printing/str.py
index 63c483c..615c2f0 100644
--- a/sympy/printing/str.py
+++ b/sympy/printing/str.py
@@ -204,6 +204,11 @@ class StrPrinter(Printer):
         else:
             return 'O(%s)'%self.stringify(expr.args, ', ', 0)
 
+    def _print_PDF(self, expr):
+        return 'PDF(%s, (%s, %s, %s))' % \
+            (self._print(expr.pdf.args[1]), self._print(expr.pdf.args[0]), \
+            self._print(expr.domain[0]), self._print(expr.domain[1]))
+
     def _print_Pi(self, expr):
         return 'pi'
 
diff --git a/sympy/statistics/distributions.py b/sympy/statistics/distributions.py
index f8ebc8e..5903400 100644
--- a/sympy/statistics/distributions.py
+++ b/sympy/statistics/distributions.py
@@ -208,3 +208,103 @@ class Uniform(ContinuousProbability):
         m = sample.mean
         d = sqrt(12*sample.variance)/2
         return Uniform(m-d, m+d)
+
+class PDF(ContinuousProbability):
+    def __init__(self, pdf, var):
+        #XXX maybe add some checking of parameters
+        if  isinstance(var, (tuple, list)):
+            self.pdf = Lambda(var[0], pdf)
+            self.domain = tuple(var[1:])
+        else:
+            self.pdf = Lambda(var, pdf)
+            self.domain = (-oo, oo)
+        self._cdf = None
+        self._mean = None
+        self._variance = None
+        self._stddev = None
+
+
+    def normalize(self):
+        norm = self.probability(*self.domain)
+        if norm != 1:
+            w = Symbol('w', real=True, dummy=True)
+            return self.__class__(self.pdf(w)/norm, (w, self.domain[0], self.domain[1]))
+            #self._cdf = Lambda(w, (self.cdf(w) - self.cdf(self.domain[0]))/norm)
+            #if self._mean is not None:
+            #    self._mean /= norm
+            #if self._variance is not None:
+            #    self._variance = (self._variance + (self._mean*norm)**2)/norm - self.mean**2
+            #if self._stddev is not None:
+            #    self._stddev = sqrt(self._variance)
+        else:
+            return self
+
+
+    def cdf(self, x):
+        x = sympify(x)
+        if self._cdf is not None:
+            return self._cdf(x)
+        else:
+            from sympy import integrate
+            w = Symbol('w', real=True, dummy=True)
+            self._cdf = integrate(self.pdf(w), w)
+            self._cdf = Lambda(w, self._cdf - self._cdf.subs(w, self.domain[0]))
+            return self._cdf(x)
+
+    def _get_mean(self):
+        if self._mean is not None:
+            return self._mean
+        else:
+            from sympy import integrate
+            w = Symbol('w', real=True, dummy=True)
+            self._mean = integrate(self.pdf(w)*w,(w,self.domain[0],self.domain[1]))
+            return self._mean
+
+    def _get_variance(self):
+        if self._variance is not None:
+            return self._variance
+        else:
+            from sympy import integrate, trim, together
+            w = Symbol('w', real=True, dummy=True)
+            self._variance = integrate(self.pdf(w)*w**2,(w,self.domain[0],self.domain[1])) - self.mean**2
+            self._variance = trim(self._variance)
+            return self._variance
+
+    def _get_stddev(self):
+        if self._stddev is not None:
+            return self._stddev
+        else:
+            self._stddev = sqrt(self.variance)
+            return self._stddev
+
+    mean = property(_get_mean)
+    variance = property(_get_variance)
+    stddev = property(_get_stddev)
+
+
+    def _random(s):
+        raise NotImplementedError
+
+    def transform(self,func,var):
+        """Return a probability distribution of random variable func(x)"""
+        #gamma = func(xi)
+
+        #w1 = Symbol('w1', real=True, dummy=True)
+        w = Symbol('w', real=True, dummy=True)
+
+        from sympy import solve
+        inverse = solve(func-w, var)
+        newPdf = S.Zero
+        funcdiff = func.diff(var)
+        #TODO check if x is in domain
+        for x in inverse:
+            # this assignment holds only for x in domain
+            # in general it would require implementing
+            # piecewise defined functions in sympy
+            newPdf += (self.pdf(var)/abs(funcdiff)).subs(var,x)
+        #from sympy import pprint
+        #pprint( newPdf)
+
+        #TODO compute new domain
+
+        return PDF(newPdf, (w, func.subs(var, self.domain[0]), func.subs(var, self.domain[1])))
-- 
2.11.4.GIT