w1 = 1.5
w2 = 1.0
u1 = 0.7
u2 = 0.6
u3 = 0.0
N = 100000
M = array([[w1 * (1 - u1), w1 * u3], [w1 * u1, w2 * (1 - u2)]]) / M[1,1]
print "M:\n",M
s = 1 - w2*(1-u2)/w1
print s, u1, exp(-u1/s), u1/s * exp(-u1/s)
M: [[ 0.0757809 0. ] [ 0.17682209 0.0673608 ]] 0.733333333333 0.7 0.384987098923 0.367487685336
max_t = 50
p = [None] * max_t
p[0] = array([1.,0.])
for t in range(1,max_t):
#p[t] = array([w1*(1-u1)*p[t-1][0] + w1*u3*p[t-1][1] , w2*u1*p[t-1][0] + w2*(1-u2)*p[t-1][1]])
p[t] = M.dot(p[t-1])
p[t] /= p[t].sum()
p = array(p)
plot(range(max_t), (p[:,0]))#/p.sum(axis=1)))
xlabel('Generations')
ylabel('A1 frequency')
ylim(0,1)
axhline(1-u1/s);
fig,ax = plt.subplots(1,1)
ax.plot(linspace(0,1), 3*(1-2*linspace(0,1)), '-k')
ax.plot(linspace(0,1), 2*(1-linspace(0,1)), '--k')
ax.set_xlabel(r'$\mu$')
ax.set_ylabel('Replacement rate')
ax.set_yticks([3,2,1])
ax.set_xticks([])
ax.set_yticklabels(['w1','w2','1'])
ax.set_xlim(0,0.5)
ax.set_ylim(1,3.25)
ax.set_title("Figure 3 left");
w2
1.0