Repository: sympy/sympy

Issue:
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]])
```

Discussion:

Patch:
diff --git a/sympy/stats/stochastic_process_types.py b/sympy/stats/stochastic_process_types.py
index 0de4d3eeaa..fdac9cb8cc 100644
--- a/sympy/stats/stochastic_process_types.py
+++ b/sympy/stats/stochastic_process_types.py
@@ -1055,13 +1055,32 @@ def communication_classes(self) -> tList[tTuple[tList[Basic], Boolean, Integer]]
         return sympify(list(zip(classes, recurrence, periods)))
 
     def fundamental_matrix(self):
+        """
+        Each entry fundamental matrix can be interpreted as
+        the expected number of times the chains is in state j
+        if it started in state i.
+
+        References
+        ==========
+
+        .. [1] https://lips.cs.princeton.edu/the-fundamental-matrix-of-a-finite-markov-chain/
+
+        """
         _, _, _, Q = self.decompose()
-        if Q == None:
-            return None
-        I = eye(Q.shape[0])
-        if Q.shape[0]*Q.shape[1] == 0 or (I - Q).det() == 0:
-            raise ValueError("Fundamental matrix doesn't exists.")
-        return ImmutableMatrix((I - Q).inv().tolist())
+
+        if Q.shape[0] > 0:  # if non-ergodic
+            I = eye(Q.shape[0])
+            if (I - Q).det() == 0:
+                raise ValueError("The fundamental matrix doesn't exist.")
+            return (I - Q).inv().as_immutable()
+        else:  # if ergodic
+            P = self.transition_probabilities
+            I = eye(P.shape[0])
+            w = self.fixed_row_vector()
+            W = Matrix([list(w) for i in range(0, P.shape[0])])
+            if (I - P + W).det() == 0:
+                raise ValueError("The fundamental matrix doesn't exist.")
+            return (I - P + W).inv().as_immutable()
 
     def absorbing_probabilities(self):
         """


Test patch:
diff --git a/sympy/stats/tests/test_stochastic_process.py b/sympy/stats/tests/test_stochastic_process.py
index 24ae6db694..50511d9dea 100644
--- a/sympy/stats/tests/test_stochastic_process.py
+++ b/sympy/stats/tests/test_stochastic_process.py
@@ -120,7 +120,7 @@ def test_DiscreteMarkovChain():
     TO3 = Matrix([[Rational(1, 4), Rational(3, 4), 0],[Rational(1, 3), Rational(1, 3), Rational(1, 3)], [0, Rational(1, 4), Rational(3, 4)]])
     Y2 = DiscreteMarkovChain('Y', trans_probs=TO2)
     Y3 = DiscreteMarkovChain('Y', trans_probs=TO3)
-    raises (ValueError, lambda: Y3.fundamental_matrix())
+    assert Y3.fundamental_matrix() == ImmutableMatrix([[176, 81, -132], [36, 141, -52], [-44, -39, 208]])/125
     assert Y2.is_absorbing_chain() == True
     assert Y3.is_absorbing_chain() == False
     assert Y2.canonical_form() == ([0, 1, 2], TO2)
@@ -141,6 +141,12 @@ def test_DiscreteMarkovChain():
     Y6 = DiscreteMarkovChain('Y', trans_probs=TO6)
     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)]])
     assert Y6.absorbing_probabilities() == ImmutableMatrix([[Rational(3, 4), Rational(1, 4)], [S.Half, S.Half], [Rational(1, 4), Rational(3, 4)]])
+    TO7 = Matrix([[Rational(1, 2), Rational(1, 4), Rational(1, 4)], [Rational(1, 2), 0, Rational(1, 2)], [Rational(1, 4), Rational(1, 4), Rational(1, 2)]])
+    Y7 = DiscreteMarkovChain('Y', trans_probs=TO7)
+    assert Y7.is_absorbing_chain() == False
+    assert Y7.fundamental_matrix() == ImmutableDenseMatrix([[Rational(86, 75), Rational(1, 25), Rational(-14, 75)],
+                                                            [Rational(2, 25), Rational(21, 25), Rational(2, 25)],
+                                                            [Rational(-14, 75), Rational(1, 25), Rational(86, 75)]])
 
     # test for zero-sized matrix functionality
     X = DiscreteMarkovChain('X', trans_probs=Matrix([[]]))