time = 20
numTrials = 1000
X = np.empty( [time, numTrials])
# T is a transition probability matrix for a directed graph on 4 vertices
T = np.array([[ 0, 0, 1, 0.5 ],[ 1/3, 0, 0, 0 ],[ 1/3, 0.5, 0, 0.5],[ 1/3, 0.5, 0, 0 ] ] )
for trial in np.arange(0, numTrials ):
X[0, trial] = np.random.choice( np.arange(0, 4) )
for k in np.arange(1, time):
currentPosition = X[ k-1, trial].astype( int )
X[ k, trial ] = np.random.choice( np.arange(0, 4) , p = T[:, currentPosition] )
finalPositions = X[ time-1, : ].astype(int)
frequency = np.bincount(finalPositions)
distribution = frequency/np.sum(frequency)
display(distribution)