lambda0,lambda1,lambda2 = 3,5,7
m = 1
c1,c2 = 10,12

# In scipy, you can find a function pdtr(c,rho) = Prob(Poisson(rho)<=c).
# We'll define dois(c,rho) = Prob(Poisson(rho)=c) ourselves from it.
import scipy.special, numpy
def dpois(c,rho): return scipy.special.pdtr(c,rho)-scipy.special.pdtr(c-1,rho)
def erlang(rho,c): return dpois(c,rho)/scipy.special.pdtr(c,rho)


def updateB((B1old,B2old)):
    B1new = erlang((lambda1+lambda0*(1-B2old))*m,c1)
    B2new = erlang((lambda2+lambda0*(1-B1new))*m,c2)
    return (B1new,B2new)

# I'll store successive estimates for (B1,B2) in a list.
# Start off with some initial guess e.g. (.5,.5)
# Repeatedly take the last guess (the last element in the list) and update it
# and attach this new estimate to the end of the list.
bvals = [ (.5,.5) ]
for i in range(10):
    bold = bvals[-1]
    bnew = updateB(bold)
    bvals.append(bnew)
print numpy.array(bvals)