LiveBirthDeathSequences

936 days ago by raazesh.sainudiin

Number of Live or non-extinct Birth-Death Sequences

2016 Raazesh Sainudiin

A simple generator of all binary birth-death or 1-0 sequences of length n starting from one individual. Here, 1 indicates a birth event and a 0 indicates a death event. We are interested in sequences whose number of 0s is strictly less than the number of 1s for any of its subsequence of length m <= n.

The number of such sequences for a given n is captured by https://oeis.org/A001405 with a(n) = binomial(n, floor(n/2)).

# All possible live or non-extinct birth-death or 1-0 sequences of length n # # n non-extinct birth-death sequences # ------------------------------------------ # 0 [[1.0]] # 1 [[1.0, 1.0]] # 2 [[1.0, 1.0, 0.0], [1.0, 1.0, 1.0]] # 3 [[1.0, 1.0, 0.0, 1.0], [1.0, 1.0, 1.0, 0.0], [1.0, 1.0, 1.0, 1.0]] 
       

The following paper gives the number of such sequences that terminate at a given cummulative sum (here 0s are replaced by -1s WLOG):

http://www.maa.org/sites/default/files/D11233._F._Bailey.pdf

$$A(n,k) = \frac{(n+1-k)(n+2)(n+3)\cdots(n+k)}{k!}$$

Here $A(n,k)$ denotes the number of $n$ $+1$'s and $k$ $-1$'s in a sequence

$$a_1,a_2,\ldots, a_{n+k}, \text{such that}, \quad \sum_{i=1}^m a_i \geq 0, \forall 1 \leq m \leq n + k $$

def A(n,k): return (n+1-k)*prod([x for x in range(n+2,n+k+1)])/factorial(k) 
       
T=5; a= [(T-k,k,A(T-k,k)) for k in range(0,(T/2)+1)]; a 
       
[(5, 0, 6), (4, 1, 4), (3, 2, 5)]
[(5, 0, 6), (4, 1, 4), (3, 2, 5)]
import numpy as np def BD(b,d): ''' from https://en.wikipedia.org/wiki/Bertrand's_ballot_theorem; see Feller vol 1 Bertrand's proof shows that the number of favourable orders in which to count the votes is (though he does not give this number explicitly)''' return binomial(b+d-1,b-1)-binomial(b+d-1,b) 
       
# not useful yet as it misses boundary conditions var ('B','D'); assume(B,'integer'); assume(D, 'integer'); assume(B>D); assume(D>=0); f= binomial(B+D-1,B-1)-binomial(B+D-1,B); f.full_simplify() 
       
(B - D)*factorial(B + D - 1)/(B*D*factorial(B - 1)*factorial(D - 1))
(B - D)*factorial(B + D - 1)/(B*D*factorial(B - 1)*factorial(D - 1))
def BD0(b,d): return (b-d)*factorial(b+d-1)/(b*d*factorial(b-1)*factorial(d-1)) 
       
# starting from (0,0) includes original birth T=10; a= [(T-k,k,BD(T-k,k)) for k in range(0,(T/2)+1)]; a 
       
[(10, 0, 1), (9, 1, 8), (8, 2, 27), (7, 3, 48), (6, 4, 42), (5, 5, 0)]
[(10, 0, 1), (9, 1, 8), (8, 2, 27), (7, 3, 48), (6, 4, 42), (5, 5, 0)]
# total number of live b-d paths and T-1 are events after original birth binomial(T-1,floor((T-1)/2)) 
       
126
126
sum([x[2] for x in a]) 
       
126
126
def EqualProdProbs(T): a = [[T-k,k,BD(T-k,k)] for k in range(0,(T/2)+1)] s = sum([x[2] for x in a]) for i in range(len(a)): a[i][2]=1.0*a[i][2]/s return a 
       
EqualProdProbs(20) 
       
[[20, 0, 0.0000108250882244690],
 [19, 1, 0.000194851588040443],
 [18, 2, 0.00164541341011929],
 [17, 3, 0.00863842040312628],
 [16, 4, 0.0314685314685315],
 [15, 5, 0.0839160839160839],
 [14, 6, 0.167832167832168],
 [13, 7, 0.251748251748252],
 [12, 8, 0.272727272727273],
 [11, 9, 0.181818181818182],
 [10, 10, 0.000000000000000]]
[[20, 0, 0.0000108250882244690],
 [19, 1, 0.000194851588040443],
 [18, 2, 0.00164541341011929],
 [17, 3, 0.00863842040312628],
 [16, 4, 0.0314685314685315],
 [15, 5, 0.0839160839160839],
 [14, 6, 0.167832167832168],
 [13, 7, 0.251748251748252],
 [12, 8, 0.272727272727273],
 [11, 9, 0.181818181818182],
 [10, 10, 0.000000000000000]]
def unEqualProdProbs(T, p): a = [[T-k,k,p^(T-k)*(1-p)^(k)*BD(T-k,k)] for k in range(0,(T/2)+1)] s = sum([x[2] for x in a]) for i in range(len(a)): a[i][2]=1.0*a[i][2]/s return a 
       
a=unEqualProdProbs(20, 1/2); a.reverse(); a 
       
[[10, 10, 0.000000000000000],
 [11, 9, 0.181818181818182],
 [12, 8, 0.272727272727273],
 [13, 7, 0.251748251748252],
 [14, 6, 0.167832167832168],
 [15, 5, 0.0839160839160839],
 [16, 4, 0.0314685314685315],
 [17, 3, 0.00863842040312628],
 [18, 2, 0.00164541341011929],
 [19, 1, 0.000194851588040443],
 [20, 0, 0.0000108250882244690]]
[[10, 10, 0.000000000000000],
 [11, 9, 0.181818181818182],
 [12, 8, 0.272727272727273],
 [13, 7, 0.251748251748252],
 [14, 6, 0.167832167832168],
 [15, 5, 0.0839160839160839],
 [16, 4, 0.0314685314685315],
 [17, 3, 0.00863842040312628],
 [18, 2, 0.00164541341011929],
 [19, 1, 0.000194851588040443],
 [20, 0, 0.0000108250882244690]]
# WLLN: for p > 1/2, # (\sum_{t=1}^T X_i)/T = the fraction of live individuals by time-step T, # converges to 2p-1 as T approaches infinity # E(X_i) = (+1)*p + (-1)*(1-p) = p-1+p = 2p-1, # where X_i is +1 with prob p and -1 with prob 1-p with IID increment assumption # V(X_i)=4p(1-p) and Chebychev inequality proves WLLN 2*p-1. 
       
0.820000000000000
0.820000000000000
p=0.91; T=100; a=unEqualProdProbs(T, p); a.reverse(); pmf=np.array([x[2] for x in a]); F=np.cumsum(pmf) P = [(x[0]-x[1])*1./T for x in a] 
       
pmf.dot(np.array(P)) 
       
0.82399512195121982
0.82399512195121982
plt+=line(zip(P,F)) plt.show() 
       
# there is also a CLT for this concentration to 2p-1 as a function of N-1 
       
var('N','t','y','r') 
       
(N, t, y, r)
(N, t, y, r)
y*(N-t)/( y*(N-t)+r*y) # prob of '+1' or infection in SIR 
       
(N - t)*y/((N - t)*y + r*y)
(N - t)*y/((N - t)*y + r*y)
E = y*(N-t)/( y*(N-t)+r*y) - (1-(y*(N-t)/( y*(N-t)+r*y))) E = E.full_simplify() E 
       
(N - r - t)/(N + r - t)
(N - r - t)/(N + r - t)
E.sum(t,1,N) 
       
-2*r*(harmonic_number(-N - r) - harmonic_number(-r)) + N
-2*r*(harmonic_number(-N - r) - harmonic_number(-r)) + N
 
       
import numpy as np def String2Int(s, deLimiter): return np.array([np.double(x) for x in s.split(deLimiter)]) def isLive(s): '''return true if all subsequences have more 1's than 0's''' n=s.size seq1Ton = np.array([np.double(1.0)/np.double(x) for x in range(1,n+1)]) c = s.cumsum() * seq1Ton d = c > 1/2 return d.prod() def nextLiveSeqs(s): t=[] for x in s: #print '----- x is ------' #print x #print '-----------' y = copy(x); y.append(np.double(0.0)); if isLive(np.array(y)): t.append(y) y = copy(x); y.append(np.double(1.0)); if isLive(np.array(y)): t.append(y) #print '----- t is ------' #print t return t def compose(f, x, n): # computes f(f(...f(x))) for i in range(n): x = f(x) # this change is local to this function call! return x def zeros2minus1s(a): aa=[] for y in a: #x = np.array([np.double(1.0) for i in range(len(y))]).cumsum() yy = copy(np.array(y)) for i in range(len(yy)): if yy[i]==0: yy[i] = np.double(-1.0) yyy = np.array(yy).cumsum() #aa.append(yyy) # append np.array aa.append([ np.double(x) for x in yyy]) return aa 
       
a=EqualProdProbs(10); a.reverse(); (np.array([x[2] for x in a])).cumsum() 
       
array([ 0.16030534,  0.50381679,  0.79007634,  0.92366412,  0.95801527, 
1.        ])
array([ 0.16030534,  0.50381679,  0.79007634,  0.92366412,  0.95801527,  1.        ])
isLive(String2Int('1 1 0 1',' ')) 
       
1
1
s=[[np.double(1.0)]] s 
       
[[1.0]]
[[1.0]]
t = nextLiveSeqs([[np.double(1.0), np.double(1.0)]]) t 
       
[[1.0, 1.0, 0.0], [1.0, 1.0, 1.0]]
[[1.0, 1.0, 0.0], [1.0, 1.0, 1.0]]
zeros2minus1s(t) 
       
[[1.0, 2.0, 1.0], [1.0, 2.0, 3.0]]
[[1.0, 2.0, 1.0], [1.0, 2.0, 3.0]]
compose(nextLiveSeqs,s,0) 
       
[[1.0]]
[[1.0]]
compose(nextLiveSeqs,s,1) 
       
[[1.0, 1.0]]
[[1.0, 1.0]]
compose(nextLiveSeqs,s,2) 
       
[[1.0, 1.0, 0.0], [1.0, 1.0, 1.0]]
[[1.0, 1.0, 0.0], [1.0, 1.0, 1.0]]
compose(nextLiveSeqs,s,3) 
       
[[1.0, 1.0, 0.0, 1.0], [1.0, 1.0, 1.0, 0.0], [1.0, 1.0, 1.0, 1.0]]
[[1.0, 1.0, 0.0, 1.0], [1.0, 1.0, 1.0, 0.0], [1.0, 1.0, 1.0, 1.0]]
a = compose(nextLiveSeqs,s,5) a 
       
[[1.0, 1.0, 0.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]]
[[1.0, 1.0, 0.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]]
zeros2minus1s(compose(nextLiveSeqs,[[np.double(1.0)]],2)) 
       
[[1.0, 2.0, 1.0], [1.0, 2.0, 3.0]]
[[1.0, 2.0, 1.0], [1.0, 2.0, 3.0]]
zeros2minus1s(compose(nextLiveSeqs,[[np.double(1.0)]],3)) 
       
[[1.0, 2.0, 1.0, 2.0], [1.0, 2.0, 3.0, 2.0], [1.0, 2.0, 3.0, 4.0]]
[[1.0, 2.0, 1.0, 2.0], [1.0, 2.0, 3.0, 2.0], [1.0, 2.0, 3.0, 4.0]]
zeros2minus1s(compose(nextLiveSeqs,[[np.double(1.0)]],4)) 
       
[[1.0, 2.0, 1.0, 2.0, 1.0],
 [1.0, 2.0, 1.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 2.0, 1.0],
 [1.0, 2.0, 3.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 5.0]]
[[1.0, 2.0, 1.0, 2.0, 1.0],
 [1.0, 2.0, 1.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 2.0, 1.0],
 [1.0, 2.0, 3.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 5.0]]
zeros2minus1s(compose(nextLiveSeqs,[[np.double(1.0)]],5)) 
       
[[1.0, 2.0, 1.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]]
[[1.0, 2.0, 1.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]]
zeros2minus1s(compose(nextLiveSeqs,[[np.double(1.0)]],6)) 
       
[[1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0],
 [1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 3.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 3.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 3.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 5.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 1.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 1.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 3.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 5.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 5.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 5.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 5.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]]
[[1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0],
 [1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 3.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 3.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 3.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 5.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 1.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 1.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 3.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 5.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 5.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 5.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 5.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]]
zeros2minus1s(compose(nextLiveSeqs,[[np.double(1.0)]],7)) 
       
[[1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]]
[[1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 2.0, 3.0, 4.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0, 2.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 5.0, 4.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 5.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 6.0],
 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]]
pic = points([(0,0)]); #y=a[0] #n=len(y)-1 #pic = points([(0,0),(len(y),len(y))]); for y in a: #pic += line([(0,0),(len(y),len(y))]) pic.show(figsize=[2,2]) x = np.array([np.double(1.0) for i in range(len(y))]).cumsum() yy = copy(np.array(y)) for i in range(len(yy)): if yy[i]==0: #print 'zer' yy[i] = np.double(-1.0) yyy = np.array(yy).cumsum() pts = zip(x,yyy) #print pts g = line(pts) g += points(pts) pic += g 
       









pic.show(figsize=[2,2]) 
       
compose(nextLiveSeqs,s,6) 
       
[[1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0],
 [1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]]
[[1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0],
 [1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0],
 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]]
[len(compose(nextLiveSeqs,s,i)) for i in range(15)] 
       
[1, 1, 2, 3, 6, 10, 20, 35, 70, 126, 252, 462, 924, 1716, 3432]
[1, 1, 2, 3, 6, 10, 20, 35, 70, 126, 252, 462, 924, 1716, 3432]
[len(compose(nextLiveSeqs,s,i)) for i in range(20)] 
       
[1,
 1,
 2,
 3,
 6,
 10,
 20,
 35,
 70,
 126,
 252,
 462,
 924,
 1716,
 3432,
 6435,
 12870,
 24310,
 48620,
 92378]
[1,
 1,
 2,
 3,
 6,
 10,
 20,
 35,
 70,
 126,
 252,
 462,
 924,
 1716,
 3432,
 6435,
 12870,
 24310,
 48620,
 92378]

Next Steps

 

  1. Set up a Hasse diagram with the correct immediate precedence rules between 2D integer lattice points that are visited by such partial Dyck-path sequences of different lengths, so one can capture the growth of the sequence corresponding to new birth and death events, one at a time.
  2. Then one needs to find the exact number of ways to reach an integer pair (n,s) for each integer s > 1, recording the number of surviving individuals and each integer n > 1 recording the number of birth or death events so far. This is a problem of putting a(n)=Binomial(n,floor(n/2)) paths into n bins labelled {(n,s): s =1,2,...,n }.
  3. Get the expression for the probabilities of each sequence from the SIR model's parametrization via COntinuous time discrete state Markov chain (marginalizing over the possible transmission trees).
  4. Evaluate these expressions and see if any nice simplifications happen with regard to Pr{(n,s)}
  5. Take weak limits on 1/n {(n,s): s =1,2,...,n } as n \to \infty and hope to get the SIR dynamical system model or an appropriate projection of (n,s).