#http://stackoverflow.com/questions/22742951/solve-an-equation-using-a-python-numerical-solver-in-numpy
from scipy.optimize import fsolve
# Define the expression whose roots we want to find
N = 1e6
p0 = 0.5
pfix = 0.65
f = lambda s,N,p0 : (1.0 - exp(-2*N*p0*s))/(1.0 - exp(-2*N*s))
# Use the numerical solver to find the roots
s_initial_guess = 1.0/N
s_solution = fsolve(lambda s: pfix - f(s,N,p0), s_initial_guess)
# Plot it
s = np.logspace(-8,-5,50)
plot(s, f(N,s,p0), '--k')
xscale('log')
xlabel("$s$")
ylabel("$p_{fix}$")
plot([s_solution],[pfix],'ok')
axvline(x=s_solution, color='k')
axhline(y=pfix, color='k')
grid(False, 'minor')
print "The selection coefficient is s = %.2g" % s_solution
print "at which the fixation probability is %.2g" % f(s_solution,N,p0)
The selection coefficient is s = 6.2e-07 at which the fixation probability is 0.65
def find_s(N,p0,pfix):
''' N - effective population size
p0 - initial frequency
pfix - fixation probability
'''
s_initial_guess = 1.0 / N
return fsolve(lambda s: pfix - f(s,N,p0), s_initial_guess)
def simulation(N,p0,s):
N = float(N)
n = N*p0
while 0 < n < N:
w = (1 + s) * n/N + (N - n)/N # mean fitness
n = n*(1+s) / w # expected next generation
n = binomial(N, n/N)
return n == N
simulation(1e6,0.5,0.01)
True
def estimate_pfix(N,p0,s,reps=100):
pfix = array([simulation(N,p0,s) for _ in range(reps)])
return pfix.mean(), pfix.std(ddof=1)/sqrt(reps)
estimate_pfix(1e6,0.5,-1e-4)
(0.0, 0.0)
s_range = logspace(-4,-2,10)
pfix = array([estimate_pfix(N,p0,s) for s in s_range])
est_s = array([find_s(N,p0,p) for p in pfix])
plot(s_range, est_s, 'ok')
[<matplotlib.lines.Line2D at 0x10f91d68>]