Repository: sympy/sympy
Pull request: Fix for fundamental matrix of markov chains
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### 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 #20189 
#### Brief description of what is fixed or changed
I have updated ``` fundamental_matrix()``` function so that we can now calculate the fundamental matrix for regular markov chains. 

#### 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 support for fundamental matrix for regular markov chains
<!-- END RELEASE NOTES -->

Problem statement:
Unable to find fundamental matrix for regular markov chains
Currently, it is not possible to calculate the fundamental matrix for regular markov chains as of now ``` fundamental_matrix()``` was only implemented for the absorbing Markov chains. 
For example
```python
>>> T = Matrix([[0.5, 0.25, 0.25], [0.5, 0, 0.5], [0.25, 0.25, 0.5]])
>>> Y = DiscreteMarkovChain("Y", [0, 1, 2], T)
>>> Y.is_regular()
True
>>> Y.fundamental_matrix()
DiscreteMarkovChain(Y, FiniteSet(0, 1, 2), Matrix([
[ 0.5, 0.25, 0.25],
[ 0.5,    0,  0.5],
[0.25, 0.25,  0.5]]))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-5-bd1859c21b86> in <module>
----> 1 Y.fundamental_matrix()

~/project/sympy(main)/sympy/sympy/stats/stochastic_process_types.py in fundamental_matrix(self)
    638         I = eye(Q.shape[0])
    639         if (I - Q).det() == 0:
--> 640             raise ValueError("Fundamental matrix doesn't exists.")
    641         return ImmutableMatrix((I - Q).inv().tolist())
    642 

ValueError: Fundamental matrix doesn't exists.
```
But correct answer should be
``` 
Matrix([
[  1.14666666666667, 0.04, -0.186666666666667],
[              0.08, 0.84, 0.0799999999999999],
[-0.186666666666667, 0.04,   1.14666666666667]])
```

Hints from issue discussion:

Commit message: added support for fundamental matrix for regular markov chains

Patch to review:
diff --git a/sympy/stats/stochastic_process_types.py b/sympy/stats/stochastic_process_types.py
index 0e84402c3964..8775b0c1e497 100644
--- a/sympy/stats/stochastic_process_types.py
+++ b/sympy/stats/stochastic_process_types.py
@@ -846,9 +846,16 @@ def fundamental_matrix(self):
         if Q == None:
             return None
         I = eye(Q.shape[0])
-        if (I - Q).det() == 0:
-            raise ValueError("Fundamental matrix doesn't exists.")
-        return ImmutableMatrix((I - Q).inv().tolist())
+        if (self.is_absorbing_chain()):
+            if (I - Q).det() == 0:
+                raise ValueError("Fundamental matrix doesn't exists.")
+            return ImmutableMatrix((I - Q).inv().tolist())
+        else:
+            w = self.fixed_row_vector()
+            W = Matrix([list(w) for i in range(0, Q.shape[0])])
+            if (I - Q + W).det() == 0:
+                raise ValueError("Fundamental matrix doesn't exists.")
+            return ImmutableMatrix((I - Q + W).inv().tolist())
 
     def absorbing_probabilities(self):
         """
diff --git a/sympy/stats/tests/test_stochastic_process.py b/sympy/stats/tests/test_stochastic_process.py
index 2d621e4901ee..20b6da78513d 100644
--- a/sympy/stats/tests/test_stochastic_process.py
+++ b/sympy/stats/tests/test_stochastic_process.py
@@ -123,7 +123,6 @@ def test_DiscreteMarkovChain():
     Y2 = DiscreteMarkovChain('Y', trans_probs=TO2)
     Y3 = DiscreteMarkovChain('Y', trans_probs=TO3)
     assert Y3._transient2absorbing() == None
-    raises (ValueError, lambda: Y3.fundamental_matrix())
     assert Y2.is_absorbing_chain() == True
     assert Y3.is_absorbing_chain() == False
     TO4 = Matrix([[Rational(1, 5), Rational(2, 5), Rational(2, 5)], [Rational(1, 10), S.Half, Rational(2, 5)], [Rational(3, 5), Rational(3, 10), Rational(1, 10)]])


Code review comments:
diff_hunk: @@ -846,9 +846,16 @@ def fundamental_matrix(self):
         if Q == None:
             return None
         I = eye(Q.shape[0])
-        if (I - Q).det() == 0:
-            raise ValueError("Fundamental matrix doesn't exists.")
-        return ImmutableMatrix((I - Q).inv().tolist())
+        if (self.is_absorbing_chain()):
+            if (I - Q).det() == 0:
+                raise ValueError("Fundamental matrix doesn't exists.")
+            return ImmutableMatrix((I - Q).inv().tolist())
path: sympy/stats/stochastic_process_types.py
text: @user1:
have you tried `(I - Q).inv().as_immutable()`?

@author:
I am not familiar with the matrix module so can you please tell what this `as_immutable` function does exactly.
I tried modified my function by using the above lines, and it gives the same output and also passed all the test

@user2:
The changes you have made may still not work for chains which are both absorbing and regular. What should be done in that case? For this patch, absorbing chain will be preferred and even if the regular chain property is true, `ValueError` will be raised if (I - Q)^-1 doesn't exist.

@user2:
The shape of W may not always be the same as I and Q. 
```python
____________________________________________________ sympy/stats/tests/test_stochastic_process.py:test_DiscreteMarkovChain _____________________________________________________
Traceback (most recent call last):
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/stats/tests/test_stochastic_process.py", line 140, in test_DiscreteMarkovChain
    assert Y6.fundamental_matrix() == ImmutableMatrix([[Rational(3, 2), S.One, S.Half], [S.One, S(2), S.One], [S.Half, S.One, Rational(3, 2)]])
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/stats/stochastic_process_types.py", line 857, in fundamental_matrix
    if (I - Q + W).det() == 0:
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/decorators.py", line 136, in binary_op_wrapper
    return f(self)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/decorators.py", line 137, in binary_op_wrapper
    return func(self, other)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/matrices/common.py", line 2779, in __radd__
    return self + other
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/core/decorators.py", line 137, in binary_op_wrapper
    return func(self, other)
  File "/home/czgdp1807ssd/sympy_project/sympy/sympy/matrices/common.py", line 2545, in __add__
    self.shape, other.shape))
sympy.matrices.common.ShapeError: Matrix size mismatch: (3, 5) + (3, 3)
```

@author:
> The changes you have made may still not work for chains which are both absorbing and regular. What should be done in that case?

This case is not possible as chains cannot be absorbing and regular at the same time as regular matrix is a subset of ergodic chains, and also by the definition of regular chains, for some power of transition matrix it only has positive elements. 

If any example is there of this type of chains then we should add a new condition to tackle it.



@author:
> The shape of W may not always be the same as I and Q.

What I thought regarding the shape of W is that as this is for regular chains so it will have no absorbing state, so the size of W is just the sum of all the transient and recurrent state of the chains which is the shape of Q.

I have explicitly formed a W matrix from a fixed row vector just by repeating. Is there any other way of using only a fixed row vector to calculate the output i.e. something like numpy broadcasting feature?

@user2:
So, I think there shouldn't be any problem in keeping both the definitions in a single function as with any chain only one will be used? Or is there any example of a chain which is both ergodic and absorbing, though I don't think so that there should be one but in case. 

@user3:
There are chains that are neither absorbing nor ergodic. 
```
0 1 0
1 0 0
1 0 0
```

But for all possible chains, the typical definition of the fundamental matrix is needed. This is used to calculate first exit probabilities, limiting transient->recurrent probabilities, and expected number of steps within transient states.

There are cases where there are no transient states (an ergodic chain is one such instance). In this case the above 3 calculations are fairly trivial because they contain 0x0 matrices in their calculations. However, I do not believe a special exception for ergodic chains is needed (especially since there are other chains that are not ergodic that also have no transient states). I think a different function is needed for this different definition. 

This new definition is applicable only for ergodic chains while the typical definition is needed for all chains. 

@user4:
> I think there shouldn't be any problem in keeping both the definitions in a single function

Is it useful to do that though? If you were going to use the matrix for something is there a use case where either definition would be suitable?

@user2:
At most, only one definition will be used for an object, so I think, there won't be much ambiguity. The use case depends on the user, they can interpret the fundamental matrix in whatever way they want depending upon the type of the chain. If the chain is neither absorbing  nor ergodic value error should be raised. It's not a good idea IMHO to change the public APIs for this task since things are very deterministic and involve no or very little ambiguity.