From 0b1cb6ea2f696218ffbce401cbfd1b306cade428 Mon Sep 17 00:00:00 2001
From: Daniel Tomlinson <dtomlinson@panaetius.co.uk>
Date: Fri, 12 Jul 2019 01:56:47 +0100
Subject: [PATCH] updated Markov scripts for Jack

---
 markov/markov.py | 109 +++++++++++++++++++++++++----------------------
 1 file changed, 58 insertions(+), 51 deletions(-)

diff --git a/markov/markov.py b/markov/markov.py
index 976bd4b..9facfde 100644
--- a/markov/markov.py
+++ b/markov/markov.py
@@ -2,6 +2,17 @@ import numpy as np
 import itertools
 import functools
 
+
+""" Set our parameters """
+# Should we seed the results?
+setSeed = False
+seedNum = 27
+
+""" Simulation parameters """
+# Should we simulate more than once?
+setSim = True
+simNum = 1
+
 """ Define our data """
 # The Statespace
 states = np.array(['L', 'w', 'W'])
@@ -12,28 +23,23 @@ transitionName = np.array([['LL', 'Lw', 'LW'],
                            ['WL', 'Ww', 'WW']])
 
 # Probabilities Matrix (transition matrix)
-transitionMatrix = np.array([[0.8, 0.15, 0.05],
-                             [0.8, 0.15, 0.05],
-                             [0.8, 0.15, 0.05]])
+# Fill in
+transitionMatrix = None
 
 # Starting state
-startingState = 'L'
+startingState = 'w'
+initial_dist = None
+
 # Steps to run
 stepTime = 2
 # End state you want to find probabilites of
-endState = 'L'
+endState = 'W'
 
+# Get P_steps
+p_steps = False
 
-""" Set our parameters """
-# Should we seed the results?
-setSeed = False
-seedNum = 27
-
-
-""" Simulation parameters """
-# Should we simulate more than once?
-setSim = True
-simNum = 10
+# Get Stationary Dist
+stat_dist = False
 
 
 # A class that implements the Markov chain to forecast the state/mood:
@@ -65,6 +71,18 @@ class markov(object):
     def setSeed(num: int):
         return np.random.seed(num)
 
+    @staticmethod
+    def p_steps(transitionMatrix, initial_dist, steps):
+        for _ in itertools.repeat(None, steps):
+            initial_dist = transitionMatrix.T.dot(initial_dist)
+        return initial_dist
+
+    @staticmethod
+    def stationary_dist(transitionMatrix, initial_dist, steps):
+        for _ in itertools.repeat(None, steps):
+            initial_dist = transitionMatrix.T.dot(initial_dist)
+        return initial_dist
+
     @functools.lru_cache(maxsize=128)
     def forecast(self):
         print(f'Start state: {self.currentState}')
@@ -91,46 +109,21 @@ class markov(object):
                     self.currentState = "W"
                     self.stateList.append("W")
             elif self.currentState == "w":
-                self.change = np.random.choice(self.transitionName[1],
-                                               replace=True,
-                                               p=transitionMatrix[1])
-                if self.change == "ww":
-                    self.prob = self.prob * 0.15
-                    self.stateList.append("w")
-                    pass
-                elif self.change == "wL":
-                    self.prob = self.prob * 0.8
-                    self.currentState = "L"
-                    self.stateList.append("L")
-                else:
-                    self.prob = self.prob * 0.05
-                    self.currentState = "W"
-                    self.stateList.append("W")
+                # Fill in
+                pass
             elif self.currentState == "W":
-                self.change = np.random.choice(self.transitionName[2],
-                                               replace=True,
-                                               p=transitionMatrix[2])
-                if self.change == "WW":
-                    self.prob = self.prob * 0.05
-                    self.stateList.append("W")
-                    pass
-                elif self.change == "WL":
-                    self.prob = self.prob * 0.8
-                    self.currentState = "L"
-                    self.stateList.append("L")
-                else:
-                    self.prob = self.prob * 0.15
-                    self.currentState = "w"
-                    self.stateList.append("w")
+                # Fill in
+                pass
             i += 1
         print(f'Possible states: {self.stateList}')
         print(f'End state after {self.steps} steps: {self.currentState}')
         print(f'Probability of all the possible sequence of states:'
-              f' {self.prob}')
+              f' {self.prob}\n')
         return self.stateList
 
 
 def main(*args, **kwargs):
+    global startingState
     try:
         simNum = kwargs['simNum']
     except KeyError:
@@ -149,11 +142,20 @@ def main(*args, **kwargs):
     count = 0
     # Simulate Multiple Times
     if setSim:
-        for _ in itertools.repeat(None, simNum):
-            markovChain = markov(states, transitionName,
-                                 transitionMatrix, startingState,
-                                 stepTime)
-            list_state.append(markovChain.forecast())
+        if initial_dist is not None:
+            startingState = np.random.choice(states, p=initial_dist)
+            for _ in itertools.repeat(None, simNum):
+                markovChain = markov(states, transitionName,
+                                     transitionMatrix, startingState,
+                                     stepTime)
+                list_state.append(markovChain.forecast())
+                startingState = np.random.choice(states, p=initial_dist)
+        else:
+            for _ in itertools.repeat(None, simNum):
+                markovChain = markov(states, transitionName,
+                                     transitionMatrix, startingState,
+                                     stepTime)
+                list_state.append(markovChain.forecast())
     else:
         for _ in range(1, 2):
             list_state.append(markov(states, transitionName,
@@ -169,6 +171,11 @@ def main(*args, **kwargs):
         simNum = 1
     print(f'\nThe probability of starting in {startingState} and finishing'
           f' in {endState} after {stepTime} steps is {(count / simNum):.2%}')
+    if p_steps:
+        print(f'P_{stepTime} is '
+              f'{markov.p_steps(transitionMatrix, initial_dist, stepTime)}')
+    if stat_dist:
+        print(f'Stat dist is {markov.stationary_dist(transitionMatrix,initial_dist, stepTime)}')
 
 
 if __name__ == '__main__':