Python code awemany vs. dgenr8 block depth probability

stickygate.py

#!/usr/bin/env python
import numpy as np
from scipy.special import binom
from scipy.integrate import quad
from scipy.stats import erlang
######################################################################
# awemany: look at binomial
def bprop(k, n, p):
    """Binomial probability density.
    Probability of drawing exactly k x-like items out of n total, where
    probability of an item having x is p.
    """
    return binom(n, k)*(p**k)*(1-p)**(n-k)

def prob_sticky_when_x_blocks(K, a):
    """ Probability that sticky flag is reached in a set of K blocks. """
    p=0.0 # total probability
    for k in range(K, 2*K):
        # in each step, add probability that k out of x are marked with
        # probability per item a (hash power fraction)
        p+=bprop(k, 2*K-1, a)
    return p
######################################################################
# dgenr8: look at times
def pdf_to_produce_K_blocks(K, a):
    """Returns PDF for time in seconds to produce K excessive blocks with hash
    power fraction a."""
    def pdf(t):
        return erlang.pdf(t, a=K, scale=600.0/a)
    return np.vectorize(pdf)

def prob_sticky_when_x_blocks_dgenr8(K, a):
    pdf_a=pdf_to_produce_K_blocks(K, a)
    pdf_b=pdf_to_produce_K_blocks(K, 1.0-a)

    # integral routine of scipy complains about instability - let's rather use
    # the ready-made CDF
    #def cdf_b(t):
    #    return quad(pdf_b, 0.0, t)[0]
    
    cdf_a=lambda t : erlang.cdf(t, a=K, scale=600.0/a)
    cdf_b=lambda t : erlang.cdf(t, a=K, scale=600.0/(1.0-a))

    # similarly, integrating until infty fails w/ scipy - assume 10*K*600/a/(1.0-a) is
    # far enough into the tail.
    return quad(lambda t : pdf_b(t) * cdf_a(t), 0.0, 10*K*600/a/(1.0-a))[0]
######################################################################

print "K q Pawemany Pdgenr8 Pawemany-Pdgenr8"
for K in range(1, 10):
    for q in [0.05, 0.10, 0.15, 0.20, 0.25, 0.30,
              0.35, 0.40, 0.45, 0.50, 0.55, 0.60]:
        Pawemany=prob_sticky_when_x_blocks(K, q)
        Pdgenr8 =prob_sticky_when_x_blocks_dgenr8(K, q)
        print "%3d %4.2f %10.8f %10.8f %g" % (
            K, q, Pawemany, Pdgenr8, Pawemany-Pdgenr8)
    print

stickygate.txt

K q Pawemany Pdgenr8 Pawemany-Pdgenr8
  1 0.05 0.05000000 0.05000000 6.93889e-18
  1 0.10 0.10000000 0.10000000 0
  1 0.15 0.15000000 0.15000000 0
  1 0.20 0.20000000 0.20000000 0
  1 0.25 0.25000000 0.25000000 0
  1 0.30 0.30000000 0.30000000 3.27516e-15
  1 0.35 0.35000000 0.35000000 3.90465e-13
  1 0.40 0.40000000 0.40000000 1.3888e-11
  1 0.45 0.45000000 0.45000000 2.23363e-10
  1 0.50 0.50000000 0.50000000 2.06115e-09
  1 0.55 0.55000000 0.54999999 1.2698e-08
  1 0.60 0.60000000 0.59999994 5.77775e-08

  2 0.05 0.00725000 0.00725000 1.73472e-18
  2 0.10 0.02800000 0.02800000 3.46945e-18
  2 0.15 0.06075000 0.06075000 6.93889e-18
  2 0.20 0.10400000 0.10400000 4.16334e-17
  2 0.25 0.15625000 0.15625000 0
  2 0.30 0.21600000 0.21600000 0
  2 0.35 0.28175000 0.28175000 -5.55112e-17
  2 0.40 0.35200000 0.35200000 0
  2 0.45 0.42525000 0.42525000 1.11022e-16
  2 0.50 0.50000000 0.50000000 1.66533e-16
  2 0.55 0.57475000 0.57475000 5.9952e-15
  2 0.60 0.64800000 0.64800000 1.14686e-13

  3 0.05 0.00115813 0.00115812 4.33681e-19
  3 0.10 0.00856000 0.00856000 3.46945e-18
  3 0.15 0.02661187 0.02661187 0
  3 0.20 0.05792000 0.05792000 4.16334e-17
  3 0.25 0.10351562 0.10351562 0
  3 0.30 0.16308000 0.16308000 -5.55112e-17
  3 0.35 0.23516937 0.23516937 0
  3 0.40 0.31744000 0.31744000 1.66533e-16
  3 0.45 0.40687313 0.40687312 1.66533e-16
  3 0.50 0.50000000 0.50000000 1.11022e-16
  3 0.55 0.59312688 0.59312688 0
  3 0.60 0.68256000 0.68256000 0

  4 0.05 0.00019358 0.00019358 -5.42101e-20
  4 0.10 0.00272800 0.00272800 0
  4 0.15 0.01210317 0.01210317 -8.67362e-18
  4 0.20 0.03334400 0.03334400 1.38778e-17
  4 0.25 0.07055664 0.07055664 0
  4 0.30 0.12603600 0.12603600 -1.11022e-16
  4 0.35 0.19984573 0.19984573 -1.11022e-16
  4 0.40 0.28979200 0.28979200 0
  4 0.45 0.39171220 0.39171220 1.11022e-16
  4 0.50 0.50000000 0.50000000 -1.11022e-16
  4 0.55 0.60828780 0.60828780 -2.22045e-16
  4 0.60 0.71020800 0.71020800 0

  5 0.05 0.00003322 0.00003322 3.52366e-19
  5 0.10 0.00089092 0.00089092 2.56956e-17
  5 0.15 0.00562866 0.00562866 8.67362e-19
  5 0.20 0.01958144 0.01958144 4.64559e-15
  5 0.25 0.04892731 0.04892731 1.38778e-17
  5 0.30 0.09880866 0.09880866 -4.16334e-17
  5 0.35 0.17171929 0.17171929 8.32667e-17
  5 0.40 0.26656768 0.26656768 1.11022e-16
  5 0.45 0.37857905 0.37857905 2.77556e-16
  5 0.50 0.50000000 0.50000000 1.11022e-16
  5 0.55 0.62142095 0.62142095 0
  5 0.60 0.73343232 0.73343232 -1.11022e-16

  6 0.05 0.00000580 0.00000580 2.24379e-18
  6 0.10 0.00029571 0.00029571 2.46168e-16
  6 0.15 0.00265686 0.00265686 1.73472e-18
  6 0.20 0.01165421 0.01165421 5.20417e-18
  6 0.25 0.03432751 0.03432751 -2.08167e-17
  6 0.30 0.07822479 0.07822479 -8.32667e-17
  6 0.35 0.14868372 0.14868372 -5.55112e-17
  6 0.40 0.24650187 0.24650187 -2.77556e-17
  6 0.45 0.36687742 0.36687742 1.66533e-16
  6 0.50 0.50000000 0.50000000 -2.22045e-16
  6 0.55 0.63312258 0.63312258 0
  6 0.60 0.75349813 0.75349813 0

  7 0.05 0.00000103 0.00000103 -1.2678e-18
  7 0.10 0.00009929 0.00009929 -2.88953e-16
  7 0.15 0.00126755 0.00126755 6.93889e-18
  7 0.20 0.00700356 0.00700356 6.07153e-18
  7 0.25 0.02429014 0.02429014 1.38778e-17
  7 0.30 0.06237521 0.06237521 -2.08167e-17
  7 0.35 0.12946823 0.12946823 1.11022e-16
  7 0.40 0.22884395 0.22884395 8.32667e-17
  7 0.45 0.35625819 0.35625819 3.88578e-16
  7 0.50 0.50000000 0.50000000 2.22045e-16
  7 0.55 0.64374181 0.64374181 0
  7 0.60 0.77115605 0.77115605 2.22045e-16

  8 0.05 0.00000018 0.00000018 7.38904e-19
  8 0.10 0.00003362 0.00003362 6.09864e-20
  8 0.15 0.00060961 0.00060961 -9.43256e-18
  8 0.20 0.00423975 0.00423975 8.67362e-18
  8 0.25 0.01729984 0.01729984 1.38778e-17
  8 0.30 0.05001254 0.05001254 -2.08167e-17
  8 0.35 0.11323113 0.11323113 6.93889e-17
  8 0.40 0.21310318 0.21310318 3.33067e-16
  8 0.45 0.34649607 0.34649607 3.33067e-16
  8 0.50 0.50000000 0.50000000 1.11022e-16
  8 0.55 0.65350393 0.65350393 3.33067e-16
  8 0.60 0.78689682 0.78689682 2.22045e-16

  9 0.05 0.00000003 0.00000003 2.16336e-13
  9 0.10 0.00001146 0.00001146 2.87991e-20
  9 0.15 0.00029503 0.00029503 2.27682e-18
  9 0.20 0.00258146 0.00258146 7.80626e-18
  9 0.25 0.01238478 0.01238478 2.42861e-17
  9 0.30 0.04027694 0.04027694 6.93889e-18
  9 0.35 0.09937886 0.09937886 1.11022e-16
  9 0.40 0.19893649 0.19893649 2.77556e-16
  9 0.45 0.33743562 0.33743562 1.66533e-15
  9 0.50 0.50000000 0.50000000 3.33067e-16
  9 0.55 0.66256438 0.66256438 3.33067e-16
  9 0.60 0.80106351 0.80106351 3.33067e-16