In [64]:
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import numpy as np
import scipy as sp
import scipy.misc
Neighbour marginal
$X_j \sim Bin(F, p_i)$, where $p_i = {1 \over d_i}$.
$\mathbb{E}[X_j] = F/d_i$
"neighbour j not getting any frogs"
$\mathbb{P}( X_j = 0 ) = (1 - {1 \over d_i})^F$.
"neighbour j getting some frogs"
$\mathbb{P}( X_j > 0 )
= 1 - (1 - {1 \over d_i})^F$.
The random variables $\{X_j\}_j$ are not independent. Hence, the events $\{X_j>0\}$ and $\{X_k>0\}$ are not independent. Using linearity of expectation:
"number of neighbours getting frogs"
$\mathbb{E}[|\{j: X_j>0\}|] = d_i[1 - (1 - {1 \over d_i})^F]$
$\mathbb{E}[X_j | X_j>0] = {\mathbb{E}[X_j] \over \mathbb{P}[X_j>0] } = {F \over d_i [1 - (1 - {1 \over d_i})^F] }$
In [58]:
F = np.linspace(1, 40)
d = 5
f = (d * ( 1 - (1 - 1/float(d))**F))/float(d)
plt.plot(F, f, lw=4, alpha=0.7)
d = 10
f = (d * ( 1 - (1 - 1/float(d))**F))/float(d)
plt.plot(F, f, lw=4, alpha=0.7)
d = 20
f = (d * ( 1 - (1 - 1/float(d))**F))/float(d)
plt.plot(F, f, lw=4, alpha=0.7)
d = 40
f = (d * ( 1 - (1 - 1/float(d))**F))/float(d)
plt.plot(F, f, lw=4, alpha=0.7)
#plt.plot(F, (F+1)/(2*F))
#plt.plot((F, F), (0, 1), 'k-', lw=4, alpha=0.5)
plt.axis([F[0], F[-1], 0, 1])
plt.grid()
plt.legend([5, 10, 20, 40], title='Out-degree', loc='best', framealpha=0.5)
plt.xlabel('# frogs')
plt.ylabel('fraction of successors that get frogs')
plt.title('Fraction of touched neighbours')
Out[58]:
Strictly positive marginal
Let $Y_j$ denote $X_j|_{X_j>0}$.
$\mathbb{P}(Y_j = k) = {\mathbb{P}(X_j = k) \over \mathbb{P}[X_j>0] } = { {F \choose k } p_i^k (1-p_i)^{F-k} \over 1 - (1-p_i)^F } , \quad 1 \leq k \leq F $
The distribution of $Y_j$ is not binomial anymore and does not seem to have a coin tossing interpretation.
Binomial approximation
We can get a decent approximation of $Y_j$ by $$Z_j \sim 1 + Bin(F-1, q_i),$$ where $$q_i = { {F \over d_i[1 - (1 - {1 \over d_i})^F] } - 1 \over F - 1 } $$
In [482]:
def plotcomparison2(d):
froglist = [20, 100, 200]
TV = []
for F in froglist:
k=np.array(range(1,F+1))
p_i = 1/float(d)
pexact = sp.misc.comb(F,k) * np.power(p_i, k) * np.power( (1-p_i), (F-k) )
pexact = pexact / (1 - np.power( (1-p_i), F ))
q_i = F / (d*(1 - np.power( (1-p_i), F ))) - 1
q_i = q_i / (F-1)
papprox = sp.misc.comb(F-1,k-1) * np.power(q_i, k-1) * np.power( (1-q_i), ((F-1)-(k-1) ) )
plt.plot(k, np.concatenate((pexact.reshape((F,1)), papprox.reshape((F,1))), axis=1), lw=3, alpha=0.7)
TV.append( sum(abs(pexact - papprox)) )
plt.axis([k[0], 20, 0, min([1, 2.5*max(pexact)]) ])
plt.grid()
#plt.legend(['Exact', 'Binom. Approx.'], title='Distribution', loc='best', framealpha=0.5)
plt.xlabel('k')
plt.ylabel('PMF')
plt.title(str(froglist) + ' frogs, ' + str(d) + ' neighbours (T.V. = ' + str( ["{0:.3f}".format(c) for c in TV] ) + ')' )
In [483]:
plotcomparison2(d = 20)
In [430]:
def approxfrogoutput(F,d):
k=np.array(range(1,F+1))
p_i = 1/float(d)
q_i = F / (d*(1 - np.power( (1-p_i), F ))) - 1
q_i = q_i / (F-1)
papprox = sp.misc.comb(F-1,k-1) * np.power(q_i, k-1) * np.power( (1-q_i), ((F-1)-(k-1) ) )
k=np.array(range(0,F+1))
papprox2 = np.concatenate(
(
[np.power( (1-p_i), F )],
np.array(1 - np.power( (1-p_i), F ))*papprox
)
)
papprox = papprox2
plt.semilogx(k/float(F), papprox, lw=3, alpha=0.7)
pconv = np.convolve(papprox, papprox)
plt.semilogx(np.array(range(2,2+len(pconv)))/float(F), pconv, lw=3, alpha=0.7)
lgnd = ['Binom. Approx.', 'Sum of 2']
for l in range(3,d+1):#range(d):
pconv = np.convolve(pconv, papprox)
if l in [d]:
plt.semilogx(np.array(range(0,len(pconv)))/float(F), pconv, lw=3, alpha=0.7)
lgnd.append('Sum of ' + str(l))
# Confidence interval
#cdf = np.cumsum(papprox)
rn = np.array(range(0,len(pconv)))
incl = rn[(np.cumsum(pconv) > 0.05) & (np.cumsum(pconv) < 0.95)]
shadex = rn[incl[0]:incl[-1]+2] - 0.5
shadex[shadex<1]=1
plt.fill_between(shadex/float(F), [0.0]*len(shadex), [1.0]*len(shadex), facecolor='red', alpha=0.1,
label='90% C.I.')
plt.axis([1e-1, 1e1, 0, min([1, 1.2*max(papprox)]) ])
plt.grid()
plt.legend(lgnd, title='Distribution', loc='best', framealpha=0.7)
plt.xlabel('k/F')
plt.ylabel('PMF')
plt.title('F = ' + str(F) + ' frogs, d = ' + str(d) + ' neighbours' )
In [484]:
approxfrogoutput(F=10, d = 20)
In [485]:
approxfrogoutput(F=20, d = 20)
In [486]:
approxfrogoutput(F=40, d = 20)
In [487]:
approxfrogoutput(F=80, d = 20)
In [488]:
approxfrogoutput(F=160, d = 20)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [291]:
def indepfrogoutput(F,d):
k=np.array(range(0,F+1))
p_i = 1/float(d)
pexact = sp.misc.comb(F,k) * np.power(p_i, k) * np.power( (1-p_i), (F-k) )
#pexact = pexact / (1 - np.power( (1-p_i), F ))
assert(abs(sum(pexact)-1)<1e-3)
plt.semilogx(k, pexact, lw=3, alpha=0.7)
pconv = np.convolve(pexact, pexact)
plt.semilogx(range(2,2+len(pconv)), pconv, lw=3, alpha=0.7)
lgnd = ['Binom. Approx.', 'Sum of 2']
for l in range(3,d+1):#range(d):
pconv = np.convolve(pconv, pexact)
if l in [d]:
plt.semilogx(range(0,len(pconv)), pconv, lw=3, alpha=0.7)
lgnd.append('Sum of ' + str(l))
#plt.axis([k[0], k[-1], 0, min([1, 1.5*max(pexact)]) ])
plt.grid()
plt.legend(lgnd, title='Distribution', loc='best', framealpha=0.5)
plt.xlabel('k')
plt.ylabel('PMF')
plt.title('Independent spawning: ' + str(F) + ' frogs, ' + str(d) + ' neighbours' )