from IPython.display import Image i = Image(filename='Coin_MC.png') i import numpy as np A = np.array([[.5,-.5],[-.5,1]]) b = np.array([1,1]) X = np.linalg.solve(A,b) print 'Expected number of steps from state 0: '+str(X[0]) from IPython.display import Image i = Image(filename='ABC_MC.png') i %matplotlib inline from __future__ import division import numpy as np # rates for leaving current state A = np.array([.1, .2, .1]) # number of Markov transitions to observe N = int(1e5) # state vector X = np.zeros(N+1) # time vector T = np.zeros(N+1) # 'probability' vector P = np.zeros(3) for n in range(N): # sample the next transition time from exp distribution tnext = np.random.exponential() / A[X[n]] T[n+1] = T[n] + tnext # move to the next Markov state if X[n] == 0: X[n+1] = 1 elif X[n] == 1: if np.random.rand() > .5: X[n+1] = 2 else: X[n+1] = 0 else: X[n+1] = 1; for i in range(3): # idxs of when Markov Chain transitioned to state i idx = np.where(X == i)[0]; idx = idx[:-1]; # 'Probabilities' (fraction of the time MC in state i) P[i] = np.sum( (T[idx+1] - T[idx]) ) / (T[-1] - T[0]); print "Fraction of time observed by each state is: " print P from scipy.integrate import ode def A(): ''' Transition probability matrix ''' return np.array([[-.1,.1,0], [.1,-.2,.1], [0,.1,-.1]]) def dPdt(t,P,A): ''' Differential equation: dP/dt = A * P ''' return np.dot(P,A()) def solveKFW(Tf): ''' Discretize and integrate the Kolmogorov Forward Equation. ''' #Define initial conditions P0 = np.array([1,0,0]) t0 = 0 dt = Tf/100 #set the integrator rk45 = ode(dPdt).set_integrator('dopri45') rk45.set_initial_value(P0,t0).set_f_params(A) #integrate and store each solution in #a numpy matrix, P and T. P = np.zeros((100,3)) T = np.zeros(100) P[0,:] = P0 T[0] = t0 idx = 0 while rk45.successful() and rk45.t < Tf: rk45.integrate(rk45.t+dt) P[idx,:] = np.array(rk45.y) T[idx] = rk45.t idx += 1 return P,T if __name__=='__main__': P,T = solveKFW(120) plot(T,P[:,0],'s-r') plot(T,P[:,1],'^-g') plot(T,P[:,2],'*-k') legend(['Inactive A','Active, B', 'Inactive, C']) xlabel('Time (m)') ylabel('Pr of state, s')