markov.py

#this illustrates a simple two state markov chain
#that begins all zeroes and then has transition probabilities
# of 1% in each slot; the curve shows the expected number of
# 1's in the sequence
# it's a crude model of the 4-state JC picture of evolution

import numpy as np
import matplotlib.pyplot as plt

N=3000

for trials in range(1):
    d=np.zeros(N)
    x=np.zeros(1000)
    for i in range(N):
        for j in range(len(x)):
            z=np.random.randint(0,100)
            if z==99:
                x[j]=1-x[j]
        d[i]=sum(x)/1000
    plt.plot(range(N),d)
beta=-.5*np.log(1-.02)
alpha=.5-.5*np.exp(-2*np.arange(N)*beta)
plt.plot(range(N),alpha)
plt.show()