Repository: sympy/sympy
Pull request: [GSoC] Added Poisson, Wiener and Gamma Process
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->
Added Poisson, Wiener and Gamma Process 
#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #19274 

#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* stats
  * Added Poisson, Wiener and Gamma Processes
<!-- END RELEASE NOTES -->

Problem statement:
[GSoC Discussion] Adding Poisson Process
Currently, we have DTMC, CTMC, and Bernoulli Process as the stochastic processes in `stats/stochatic_process_types.py`. Therefore, I plan to add more such processes.
We can divide(for API discussion) the stochastic processes into two categories:

**Case 1**  The processes which follow certain distribution at any time t:
(i) Poisson Process -> Poisson Distribution
(ii) Gamma Process -> Gamma Distribution
(iiI) Wiener Process -> Normal Distribution
(iv) Bernoulli Process -> Bernoulli Distribution(0/1) (Implemented)
(v) Random Walks -> Bernoulli Distribution(-1/1)

**Case 2** The processes which have no distribution:
(i) Birth Death process (should be inherited from CTMC)

So, for **Case 1** we can create a generalized framework to calculate( Expectation and Probability), in order to avoid the repetition of code in each of the classes. As, in Bernoulli Process, we will convert the Random Indexed RV to simple RV with the same name and distribution (As we already have ready-made functions for simple RV for calculating E and P) and return the E and P.

For **Case 2**, we already have P and E in CTMC and requires testing of its generalization.

ping @czgdp1807 @Upabjojr @jmig5776

Hints from issue discussion:
I'd suggest to start with designing APIs with the help of some use cases(like picking up examples from textbooks, wikipedia, research papers). In addition, some processes can be built on top of Markov chains.
For, the examples of using the processes in the **Case 1**, I have found some real-life example problems related to them [here](https://www.probabilitycourse.com/chapter11/11_1_2_basic_concepts_of_the_poisson_process.php)

Suppose if we take example 1 from the given link, then we could possibly solve it as:
```python
>>> from sympy import Eq, S
>>> from sympy.stats import PoissonProcess , P
>>> X = PoissonProcess ('X', lamda=10)
>>> P(Eq(X(S(1)/3), 2))
xxxxx         
```
In the same way, we could solve the examples of other processes too.

Therefore, the API will be the same as the API of simple RV/ simple distribution classes as in `*_types.py`, i.e, at the initialization, we need to provide the parameters that are necessary to define the distribution of the Process. The only internal difference here would be on the developer side to take care of the distribution as they vary with the time. But at the user end, we will have the same API as we have for simple RV.

For the implementation details:
https://github.com/sympy/sympy/blob/585377cc1f8876fec65968dacffd284a4572b3c9/sympy/stats/stochastic_process.py#L49-L63
As we point to the individual classes to handle the expectation and probability of the indexed RV, we need to define these methods in the respective process classes. So, to calculate E and P of expression of indexed RV using the available functions in stats, we need to convert them to simple RV with the same name and distribution and indexed RV.
So, we will globally define two functions say, `_expectation_stoch` and `_probablility_stoc`, 
(i) which will take the expression of the indexed RV, 
(ii) covert it to simple RV expression based on given process.
(iii) calculate E/P and return to the calling method of respective process classes.

Therefore all the processes in **Case 1** will commonly share two functions, which will avoid us from code repetition. 
So, it will be like:
```python

class ProcessName(...):
    def __init__(...):
        pass

    @property
    def symbol(self):
        pass

    @property
    def state_space(self):
        pass

    def distribution(self):
        pass

    def expectation(...):
        return _expectation_stoch(...)

    def probability(...):
        return _probability_stoch(...)
```
How will we solve the following part of Example 11.2 in the same link,
```
Given that the third arrival occurred at time t=2, find the probability that the fourth arrival occurs after t=4.
```

> `Given that the third arrival occurred at time t=2, find the probability that the fourth arrival occurs after t=4.`

Consider, PP as Poisson process with `lamda=2`
X1 = PP(1)
X2 = PP(2)
X05 = PP(0.5)
We can use as :
```python

>>> X05 = Poisson('X05', 2*0.5)
>>> X2 = Poisson('X2', 2*2)
>>> P(Eq(X05, 1))
exp(-1)  # Answer of 'a' there will be 1st arrival at time 0.5 secs
>>> P(Eq(X2, 0))
exp(-4)  # Answer of 'c' as there will be no arrival between t=2 and t=4, so using memory less property
```
They can be done in the above way, if we want to handle them directly as `P(Eq(X4, 4), Eq(X2, 3))`
then we need to add some additional code which will handle the memoryless property( can't be directly converted to simple RV as it will treat both of them independent and will generate wrong answer.)

Or, we can directly provide a computation of Arrival times by adding the distribution which the Arrival times follow for a particular class, and it will work fine as they are mutually independent(so they can be directly used as converting to simple RV).
What about?
```python
>>> N = PoissonProcess('N', 2)
>>> P(Eq(N(4), 4), Eq(N(2), 3))
```
Converting them to simple RV, will consider them as independent and produce an incorrect result. So, we need to add some code to apply the memory less property by which it gets converted to
```python
>>> P(Eq(N(2), 0)
```
> Converting them to simple RV

Where are they converted to simple RVs? I think we are just discussing how the interface should like. `N(t)` should be handled by `probability` to forward the call to the method of `PoissonProcess` which will solve the given query. 
I think that `probability` methods of each sub-class of stochastic processes should be able to handle the queries. If there is anything common down the hierarchy in these methods then that should be factored out and placed into the parent class.

P.S. Class design isn't a concern as of now. We should just make the interface as good as possible and then discuss about the details of implementation, that too in a PR which you will make to implement that interface. I couldn't get what's going on in https://github.com/sympy/sympy/issues/19274#issuecomment-626131730 or how the problem is solved. All the RVs and stuff should be kept hidden from the user, they should be able to make any query using `P`, and an object of `PoissonProcess` and may be if needed, some utility objects(like Matrices, assumptions). 
> Where are they converted to simple RVs? I think we are just discussing how the interface should like.

I understand I am trying to explain how we will use the available E and P if possible to use in our code. When the call will go to the method of the given process, I think this should happen,
(i) Using the given condition, we will apply memoryless property, and convert `P(expr, condition)` to `P(expr_new)` which satisfies the given conditions. We will then solve P(expr_new).

For, solving `P(expr_new)` we need to decide whether to implement new method or use the existing ones, as we had used while adding Bernoulli Process
> We should just make the interface as good as possible and then discuss about the details of implementation

The interface will allow handling the queries in the same way, as in the current implementation,
```python
>>> N = PoissonProcess('N', 2)
>>> P(Eq(N(4), 4), Eq(N(2), 3))
# This will produce required answer of 'c'
>>> P(Eq(N(0.5), 1)
# This will produce required answer of 'a'
```

Commit message: add poisson process

Patch to review:
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index cd0014de6914..7e042340b312 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -742,6 +742,10 @@ def expectation(expr, condition=None, numsamples=None, evaluate=True, **kwargs):
         evalf = kwargs.get('evalf', True)
         return sampling_E(expr, condition, numsamples=numsamples, evalf=evalf)
 
+    from sympy.stats.stochastic_process_types import StochasticProcess
+    if expr.has(StochasticProcess):
+        process = list(expr.atoms(StochasticProcess))[0]
+        return process.expectation(expr, condition, evaluate, **kwargs)
     if expr.has(RandomIndexedSymbol):
         return pspace(expr).compute_expectation(expr, condition, evaluate, **kwargs)
 
@@ -801,7 +805,10 @@ def probability(condition, given_condition=None, numsamples=None,
 
     condition = sympify(condition)
     given_condition = sympify(given_condition)
-
+    from sympy.stats.stochastic_process_types import StochasticProcess
+    if condition.has(StochasticProcess):
+        process = list(condition.atoms(StochasticProcess))[0]
+        return process.probability(condition, given_condition, evaluate, **kwargs)
     if condition.has(RandomIndexedSymbol):
         return pspace(condition).probability(condition, given_condition, evaluate, **kwargs)
 
diff --git a/sympy/stats/stochastic_process_types.py b/sympy/stats/stochastic_process_types.py
index 9380a06b92e0..fe5510252857 100644
--- a/sympy/stats/stochastic_process_types.py
+++ b/sympy/stats/stochastic_process_types.py
@@ -2,8 +2,8 @@
 
 from sympy import (Matrix, MatrixSymbol, S, Indexed, Basic,
                    Set, And, Eq, FiniteSet, ImmutableMatrix,
-                   Lambda, Mul, Dummy, IndexedBase, Add,
-                   linsolve, eye, Or, Not, Intersection,
+                   Lambda, Mul, Dummy, IndexedBase, Add, Interval, oo,
+                   linsolve, eye, Or, Not, Intersection, factorial,
                    Union, Expr, Function, exp, cacheit,
                    Ge, Piecewise, Symbol, NonSquareMatrixError)
 from sympy.core.relational import Relational
@@ -15,6 +15,7 @@
 from sympy.stats.stochastic_process import StochasticPSpace
 from sympy.stats.symbolic_probability import Probability, Expectation
 from sympy.stats.frv_types import Bernoulli, BernoulliDistribution
+from sympy.stats.drv_types import Poisson, PoissonDistribution
 from sympy.core.sympify import _sympify
 
 __all__ = [
@@ -179,7 +180,7 @@ def joint_distribution(self, *args):
                 raise ValueError("Expected a RandomIndexedSymbol or "
                                 "key not  %s"%(type(arg)))
 
-        if args[0].pspace.distribution == None: # checks if there is any distribution available
+        if not hasattr(args[0].pspace.process, 'distribution'): # checks if there is any distribution available
             return JointDistribution(*args)
 
         pdf = Lambda(tuple(args),
@@ -225,7 +226,8 @@ def __call__(self, time):
         if time not in self.index_set:
             raise IndexError("%s is not in the index set of %s"%(time, self.symbol))
         func_obj = Function(self.symbol)(time)
-        pspace_obj = StochasticPSpace(self.symbol, self)
+        distribution = getattr(self, 'distribution', None)
+        pspace_obj = StochasticPSpace(self.symbol, self, distribution)
         return RandomIndexedSymbol(func_obj, pspace_obj)
 
 class TransitionMatrixOf(Boolean):
@@ -892,11 +894,78 @@ def state_space(self):
     def distribution(self):
         return BernoulliDistribution(self.p)
 
-    def _rvindexed_subs(self, expr, condition=None):
+    def expectation(self, expr, condition=None, evaluate=True, **kwargs):
+        """
+        Computes expectation.
+
+        Parameters
+        ==========
+
+        expr: RandomIndexedSymbol, Relational, Logic
+            Condition for which expectation has to be computed. Must
+            contain a RandomIndexedSymbol of the process.
+        condition: Relational, Logic
+            The given conditions under which computations should be done.
+
+        Returns
+        =======
+
+        Expectation of the RandomIndexedSymbol.
+
+        """
+
+        return _expectation_stoch('Bernoulli', expr, condition, evaluate, **kwargs)
+
+    def probability(self, condition, given_condition=None, evaluate=True, **kwargs):
+        """
+        Computes probability.
+
+        Parameters
+        ==========
+
+        condition: Relational
+                Condition for which probability has to be computed. Must
+                contain a RandomIndexedSymbol of the process.
+        given_condition: Relational/And
+                The given conditions under which computations should be done.
+
+        Returns
+        =======
+
+        Probability of the condition.
+
+        """
+
+        return _probability_stoch('Bernoulli', condition, given_condition, evaluate, **kwargs)
+
+    def density(self, x):
+        return Piecewise((self.p, Eq(x, self.success)),
+                         (1 - self.p, Eq(x, self.failure)),
+                         (S.Zero, True))
+
+
+dict_rv_subs = {'Poisson': lambda rv: Poisson(rv.name, lamda=rv.pspace.process.lamda*rv.key),
+                'Bernoulli': lambda rv: Bernoulli(rv.name, p=rv.pspace.process.p,
+                            succ=rv.pspace.process.success, fail=rv.pspace.process.failure)
+               }
+
+
+def _rvindexed_subs(process, expr, condition=None):
         """
         Substitutes the RandomIndexedSymbol with the RandomSymbol with
         same name, distribution and probability as RandomIndexedSymbol.
 
+        Parameters
+        ==========
+
+        process : str
+           Key value of dict_rv_subs RandomVariable used to substitute
+        expr: RandomIndexedSymbol, Relational, Logic
+            Condition for which expectation has to be computed. Must
+            contain a RandomIndexedSymbol of the process.
+        condition: Relational, Logic
+            The given conditions under which computations should be done.
+
         """
 
         rvs_expr = random_symbols(expr)
@@ -904,29 +973,30 @@ def _rvindexed_subs(self, expr, condition=None):
             swapdict_expr = {}
             for rv in rvs_expr:
                 if isinstance(rv, RandomIndexedSymbol):
-                    newrv = Bernoulli(rv.name, p=rv.pspace.process.p,
-                                      succ=self.success, fail=self.failure)
+                    newrv = dict_rv_subs[process](rv)
                     swapdict_expr[rv] = newrv
             expr = expr.subs(swapdict_expr)
         rvs_cond = random_symbols(condition)
         if len(rvs_cond)!=0:
             swapdict_cond = {}
-            if condition is not None:
-                for rv in rvs_cond:
-                    if isinstance(rv, RandomIndexedSymbol):
-                        newrv = Bernoulli(rv.name, p=rv.pspace.process.p,
-                                          succ=self.success, fail=self.failure)
-                        swapdict_cond[rv] = newrv
-                condition = condition.subs(swapdict_cond)
+            for rv in rvs_cond:
+                if isinstance(rv, RandomIndexedSymbol):
+                    newrv = dict_rv_subs[process](rv)
+                    swapdict_cond[rv] = newrv
+            condition = condition.subs(swapdict_cond)
         return expr, condition
 
-    def expectation(self, expr, condition=None, evaluate=True, **kwargs):
+
+def _expectation_stoch(process, expr, condition=None, evaluate=True, **kwargs):
         """
         Computes expectation.
 
         Parameters
         ==========
 
+
+        process : str
+           Key value of dict_rv_subs RandomVariable used to substitute
         expr: RandomIndexedSymbol, Relational, Logic
             Condition for which expectation has to be computed. Must
             contain a RandomIndexedSymbol of the process.
@@ -939,7 +1009,7 @@ def expectation(self, expr, condition=None, evaluate=True, **kwargs):
         Expectation of the RandomIndexedSymbol.
 
         """
-        new_expr, new_condition = self._rvindexed_subs(expr, condition)
+        new_expr, new_condition = _rvindexed_subs(process, expr, condition)
 
         new_pspace = pspace(new_expr)
         if new_condition is not None:
@@ -951,14 +1021,15 @@ def expectation(self, expr, condition=None, evaluate=True, **kwargs):
         return new_pspace.compute_expectation(
                 new_expr, evaluate=evaluate, **kwargs)
 
-
-    def probability(self, condition, given_condition=None, evaluate=True, **kwargs):
+def _probability_stoch(process, condition, given_condition=None, evaluate=True, **kwargs):
         """
         Computes probability.
 
         Parameters
         ==========
 
+        process : str
+           Key value of dict_rv_subs RandomVariable used to substitute
         condition: Relational
                 Condition for which probability has to be computed. Must
                 contain a RandomIndexedSymbol of the process.
@@ -971,17 +1042,17 @@ def probability(self, condition, given_condition=None, evaluate=True, **kwargs):
         Probability of the condition.
 
         """
-        new_condition, new_givencondition = self._rvindexed_subs(condition, given_condition)
+        new_condition, new_givencondition = _rvindexed_subs(process, condition, given_condition)
 
         if isinstance(new_givencondition, RandomSymbol):
             condrv = random_symbols(new_condition)
             if len(condrv) == 1 and condrv[0] == new_givencondition:
-                return BernoulliDistribution(self.probability(new_condition), 0, 1)
+                return BernoulliDistribution(_probability_stoch(process, new_condition), 0, 1)
 
             if any([dependent(rv, new_givencondition) for rv in condrv]):
                 return Probability(new_condition, new_givencondition)
             else:
-                return self.probability(new_condition)
+                return _probability_stoch(process, new_condition)
 
         if new_givencondition is not None and \
                 not isinstance(new_givencondition, (Relational, Boolean)):
@@ -998,10 +1069,58 @@ def probability(self, condition, given_condition=None, evaluate=True, **kwargs):
                     % (new_condition))
         if new_givencondition is not None:  # If there is a condition
         # Recompute on new conditional expr
-            return self.probability(given(new_condition, new_givencondition, **kwargs), **kwargs)
-        return pspace(new_condition).probability(new_condition, **kwargs)
+            return _probability_stoch(process, given(new_condition, new_givencondition, **kwargs), **kwargs)
+        result = pspace(new_condition).probability(new_condition, **kwargs)
+        if evaluate and hasattr(result, 'doit'):
+            return result.doit()
+        else:
+            return result
+
+
+class PoissonProcess(ContinuousTimeStochasticProcess):
+    index_set = _set_converter(Interval(0, oo))
+
+    def __new__(cls, sym, lamda):
+        sym = _symbol_converter(sym)
+        lamda = _sympify(lamda)
+        return Basic.__new__(cls, sym, lamda)
+
+    @property
+    def state_space(self):
+        return _set_converter(Interval(0, oo))
+
+    @property
+    def symbol(self):
+        return self.args[0]
+
+    @property
+    def lamda(self):
+        return self.args[1]
+
+
+    def distribution(self, rv):
+        return PoissonDistribution(self.lamda*rv.key)
 
     def density(self, x):
-        return Piecewise((self.p, Eq(x, self.success)),
-                         (1 - self.p, Eq(x, self.failure)),
-                         (S.Zero, True))
+        return (self.lamda*x.key)**x / factorial(x) * exp(-(self.lamda*x.key))
+
+
+    def expectation(self, expr, condition=None, evaluate=True, **kwargs):
+        if given_condition is not None: # use of memory less property
+            if not isinstance(given_condition, Interval):
+                raise TypeError("Expecting an Interval as a given condition")
+            if not condition.has(PoissonProcess):
+                raise TypeError("expr should contain an instance of PoissonProcess")
+            lhs = condition.args[0](given_condition._sup - given_condition._inf)
+            condition = condition.__class__(lhs, condition.args[1])
+        return _expectation_stoch('Poisson', expr, condition=condition, evaluate=True, **kwargs)
+
+    def probability(self, condition, given_condition=None, evaluate=True, **kwargs):
+        if given_condition is not None: # use of memory less property
+            if not isinstance(given_condition, Interval):
+                raise TypeError("Expecting an Interval as a given condition")
+            if not condition.has(PoissonProcess):
+                raise TypeError("Condition should contain an instance of PoissonProcess")
+            lhs = condition.args[0](given_condition._sup - given_condition._inf)
+            condition = condition.__class__(lhs, condition.args[1])
+        return _probability_stoch('Poisson', condition, evaluate=True, **kwargs)


Code review comments:
diff_hunk: @@ -998,10 +1069,58 @@ def probability(self, condition, given_condition=None, evaluate=True, **kwargs):
                     % (new_condition))
         if new_givencondition is not None:  # If there is a condition
         # Recompute on new conditional expr
-            return self.probability(given(new_condition, new_givencondition, **kwargs), **kwargs)
-        return pspace(new_condition).probability(new_condition, **kwargs)
+            return _probability_stoch(process, given(new_condition, new_givencondition, **kwargs), **kwargs)
+        result = pspace(new_condition).probability(new_condition, **kwargs)
+        if evaluate and hasattr(result, 'doit'):
+            return result.doit()
+        else:
+            return result
+
+
+class PoissonProcess(ContinuousTimeStochasticProcess):
path: sympy/stats/stochastic_process_types.py
text: @user1:
Include it in the `sympy.stats` namespace.