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 ipywidgets import StaticInteract, RangeWidget
from ipywidgets import __version__ as ipywidgets_ver
print "ipywidgets:", ipywidgets_ver
scipy: 0.12.0 ipywidgets: 0.0.1
from mpltools import style
style.use('ggplot')
U = 0.1
s = 0.01
N = 10000
g = 25 # max_mutations
poi = poisson(U)
M = zeros((g,g))
for i in range(g-1):
for j in range(i+1):
m = poi.pmf(i-j)
M[i,j] = m
for j in range(g):
M[g-1,j] = 1-M[:,j].sum()
assert (M.sum(axis=0) == 1).all()
print M[:5,:5]
[[ 9.048e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00] [ 9.048e-02 9.048e-01 0.000e+00 0.000e+00 0.000e+00] [ 4.524e-03 9.048e-02 9.048e-01 0.000e+00 0.000e+00] [ 1.508e-04 4.524e-03 9.048e-02 9.048e-01 0.000e+00] [ 3.770e-06 1.508e-04 4.524e-03 9.048e-02 9.048e-01]]
S = diag([(1-s)**i for i in range(g)])
print S[:5,:5]
[[ 1. 0. 0. 0. 0. ] [ 0. 0.99 0. 0. 0. ] [ 0. 0. 0.98 0. 0. ] [ 0. 0. 0. 0.97 0. ] [ 0. 0. 0. 0. 0.961]]
def find_first_nonzero(x):
for i in range(g):
if x[i]>0:
return i
Start with mutation-free population:
x = zeros(g)
x[0] = 1
def ratchet_click(x):
offset = find_first_nonzero(x)
for t in range(10000):
x = M.dot(x)
x = S.dot(x)
x = x/x.sum()
assert allclose(x.sum(),1)
if N > 0:
n = multinomial(N, x)
assert len(n) == g
x = n/float(N)
if x[offset] == 0:
#print t
break
mx = arange(g).dot(x)
vx = (arange(g)**2).dot(x) - mx**2
#print "mean:", mx,U/s
#print "var:", vx,U/s
#print "mean fintess:", x.dot(S.diagonal()), exp(-U)
return x
from scipy.stats import gamma
def gamma_est(x,N,g):
counts = repeat(range(g), (x*N).astype(int))
alpha,loc,beta = gamma.fit(counts)
rv = gamma(alpha,loc=loc,scale=beta)
return rv.pdf(range(g))
def plot_ratchet(clicks):
x = zeros(g)
x[0] = 1
for i in range(clicks):
x = ratchet_click(x)
fig,ax = plt.subplots(figsize=(6, 4),
subplot_kw={'axisbg':'#EEEEEE',
'axisbelow':True})
ax.bar(range(g), x)
offset = find_first_nonzero(x)
pois = poisson(U/s).pmf(range(g))
pois2 = np.roll(pois,offset)
pois2[:offset] = 0
est = gamma_est(x,N,g)
ax.plot(range(g), pois, 'r')
ax.plot(range(g), pois2, 'g')
ax.plot(range(g), est, 'k')
ax.set_xlabel("# mutations")
ax.set_ylabel("frequency")
return fig
StaticInteract(plot_ratchet, clicks=RangeWidget(1,20,1))
c:\python27\lib\site-packages\scipy\optimize\optimize.py:438: RuntimeWarning: invalid value encountered in subtract and numpy.max(numpy.abs(fsim[0] - fsim[1:])) <= ftol): c:\python27\lib\site-packages\scipy\optimize\optimize.py:438: RuntimeWarning: invalid value encountered in absolute and numpy.max(numpy.abs(fsim[0] - fsim[1:])) <= ftol):
x = zeros(g)
x[0] = 1
N = 0
for t in range(10000):
x = M.dot(x)
x = S.dot(x)
x = x/x.sum()
assert allclose(x.sum(),1)
if N > 0:
n = multinomial(N, x)
assert len(n) == g
x = n/float(N)
fig,ax = plt.subplots(figsize=(6, 4),
subplot_kw={'axisbg':'#EEEEEE',
'axisbelow':True})
ax.bar(range(g), x, align='center')
pois = poisson(U/s).pmf(range(g))
est = gamma_est(x,1000000,g)
ax.plot(range(g), pois, 'r')
ax.plot(range(g), est, 'k')
ax.set_xlabel("# mutations")
ax.set_ylabel("frequency")
ax.set_title("MSB distribution");
Us = logspace(-2,-1)
ss = logspace(-2,-1)
n0 = array([[exp(U/s) for U in Us] for s in ss])
lvls = logspace(0,ceil(log10(n0.max())),500)
from matplotlib.colors import LogNorm
from matplotlib import ticker
figsize(8,6)
contourf(Us,ss,n0,
levels=lvls,
locator=ticker.LogLocator(),
cmap = "hot")
cb = colorbar()#format="%.2g")
cb.set_ticks(logspace(0,ceil(log10(n0.max())),1+ceil(log10(n0.max()))))
cb.ax.tick_params(labelsize=16)
cb.set_label("min N", fontsize=20)
xlabel('U', fontsize=20)
ylabel('s', fontsize=20)
xscale('log')
yscale('log')
title("min N for Poisson dist")
<matplotlib.text.Text at 0x2040feb8>
ss = logspace(-2.5,-0.5)
n0 = array([exp(U/s) for s in ss])
plot(ss,n0,'k')
xscale('log')
yscale('log')
xlabel('s', fontsize=20)
ylabel('min N', fontsize=20)
#grid(False,which='minor')
title('U=%.2g'%U)
xlim((ss.min(),ss.max()))
axvline(x=U/10, color='gray', ls='dashed')
axhline(exp(10), color='gray', ls='dashed')
text(U/9.5, exp(10.1), "%d" % exp(10));