from numpy import __version__ as np_ver
print "Numpy:", np_ver
from scipy.stats import poisson, nbinom
from scipy import __version__ as scipy_ver
print "Scipy:", scipy_ver
np.set_printoptions(linewidth=80)
np.set_printoptions(precision=3)
from IPython import __version__ as ipy_ver
print "IPython:", ipy_ver
from matplotlib import __version__ as mpl_ver
print "Matplotlib:", mpl_ver
Numpy: 1.7.1 Scipy: 0.11.0 IPython: 2.0.0-b1 Matplotlib: 1.3.1
from mpltools import style
style.use('ggplot')
--------------------------------------------------------------------------- ImportError Traceback (most recent call last) <ipython-input-8-97acd0c1ecb1> in <module>() ----> 1 from mpltools import style 2 style.use('ggplot') ImportError: No module named mpltools
Model parameters:
s = 0.1
N = 1e7
u1 = 0.003/1000
u2 = u1
x = 0.5
We need to define the drift coefficient $b(x) = E[\frac{dx}{dt}]$ and the diffusion coefficient $a(x)=E[\frac{d^2f}{dt^2}]$.
Here, $x$ is the frequency of the beneficial allele A.
The diffusion coefficient is $$ a(x) = x(1-x) $$
and the drift coefficient is $$ b(x) = Nsx(1-x)+N \mu_1(1-x)-N \mu_2 x $$
where the fitness of A is 1 and of a is 1-s, $\mu_1$ is the mutation rate from a to A, $\mu_2$ is the mutation rate from A to a, and N is the population size.
We follow ([Durret 2008(http://www.math.duke.edu/~rtd/Gbook/Gbook.html)), also see the careful derivation in this notebook.
def b(x):
return N * s * x * (1-x) + N * u1 * (1-x) - N * u2 * x
print "b",b(0.5)
def a(x):
return x * (1-x)
print "a",a(0.5)
def psi(x):
return exp(-2.0*N*s*x) * pow(x,-2.0*N*u1) * pow(1-x,-2.0*N*u2)
print "psi",psi(0.5)
from scipy.integrate import quad
def phi(x):
return quad(psi,0.5,x)[0]
print "phi",phi(0.5)
# probability to start at x and get to y before z
def P(x,y,z):
return (phi(z)-phi(x))/(phi(z)-phi(y))
print "P",P(0.5,0.6,0.4)
def m(x):
return 1.0 / (a(x) * psi(x))
# Green function
def G(x,y,u,v):
if y <= x:
return 2 * m(y) * (phi(v)-phi(x)) * (phi(y)-phi(u)) / (phi(v)-phi(u))
else:
return 2 * m(y) * (phi(x)-phi(u)) * (phi(v)-phi(y)) / (phi(v)-phi(u))
print "G",G(0.5,0.6,0.01,0.99)
def T(x,u,v):
def f(y):
return G(x,y,u,v)
return quad(f, u, v)[0]
b 250000.0 a 0.25 psi 0.0 phi 0.0 P
--------------------------------------------------------------------------- ZeroDivisionError Traceback (most recent call last) <ipython-input-12-bd07342cb9dc> in <module>() 15 def P(x,y,z): 16 return (phi(z)-phi(x))/(phi(z)-phi(y)) ---> 17 print "P",P(0.5,0.6,0.4) 18 def m(x): 19 return 1.0 / (a(x) * psi(x)) <ipython-input-12-bd07342cb9dc> in P(x, y, z) 14 # probability to start at x and get to y before z 15 def P(x,y,z): ---> 16 return (phi(z)-phi(x))/(phi(z)-phi(y)) 17 print "P",P(0.5,0.6,0.4) 18 def m(x): ZeroDivisionError: float division by zero
print "T",T(1./N,0,1)
T 3.28088687905e+27