Formula from Kimura and Ohta 1969
from scipy.integrate import romb,trapz, quad
integral = quad
def integrand1(s,N,x):
if x == 1:
return 2*N*s * (1 - exp(-2*N*s))
return (1 - exp(-2*N*s*x) - exp(-(1-x)) - exp(-2*N*s))/(x*(1-x))
def integrand2(s,N,x):
if x == 0:
return 0
return (exp(2*N*s*x) - 1) * (1 - exp(-2*N*s*x))/(x*(1-x))
def T_kimura(s,N,x):
J1 = 2/(s*(1 - exp(-2*N*s))) * integral(lambda t: integrand1(s,N,t), x, 1)[0]
u = (1-exp(-2*N*s*x))/(1-exp(-2*N*s))
J2 = 2/(s*(1-exp(-2*N*s))) * integral(lambda t: integrand2(s,N,t), 0, x)[0]
return J1 + ((1-u)/u) * J2
T_kimura(0.01, 1e5, 1e-5)