def model(args):
    '\n    %prog model erate\n\n    Model kmer distribution given error rate. See derivation in FIONA paper:\n    <http://bioinformatics.oxfordjournals.org/content/30/17/i356.full>\n    '
    from scipy.stats import binom, poisson
    p = OptionParser(model.__doc__)
    p.add_option('-k', default=23, type='int', help='Kmer size')
    p.add_option('--cov', default=50, type='int', help='Expected coverage')
    (opts, args) = p.parse_args(args)
    if (len(args) != 1):
        sys.exit((not p.print_help()))
    (erate,) = args
    erate = float(erate)
    cov = opts.cov
    k = opts.k
    xy = []
    for c in range(0, ((cov * 2) + 1)):
        Prob_Yk = 0
        for i in range((k + 1)):
            pi_i = binom.pmf(i, k, erate)
            mu_i = ((cov * ((erate / 3) ** i)) * ((1 - erate) ** (k - i)))
            Prob_Yk_i = poisson.pmf(c, mu_i)
            Prob_Yk += (pi_i * Prob_Yk_i)
        xy.append((c, Prob_Yk))
    (x, y) = zip(*xy)
    asciiplot(x, y, title='Model')
%prog model erate

Model kmer distribution given error rate. See derivation in FIONA paper:
<http://bioinformatics.oxfordjournals.org/content/30/17/i356.full>
def model(args):
    '\n    %prog model erate\n\n    Model kmer distribution given error rate. See derivation in FIONA paper:\n    <http://bioinformatics.oxfordjournals.org/content/30/17/i356.full>\n    '
    from scipy.stats import binom, poisson
    p = OptionParser(model.__doc__)
    p.add_option('-k', default=23, type='int', help='Kmer size')
    p.add_option('--cov', default=50, type='int', help='Expected coverage')
    (opts, args) = p.parse_args(args)
    if (len(args) != 1):
        sys.exit((not p.print_help()))
    (erate,) = args
    erate = float(erate)
    cov = opts.cov
    k = opts.k
    xy = []
    for c in range(0, ((cov * 2) + 1)):
        Prob_Yk = 0
        for i in range((k + 1)):
            pi_i = binom.pmf(i, k, erate)
            mu_i = ((cov * ((erate / 3) ** i)) * ((1 - erate) ** (k - i)))
            Prob_Yk_i = poisson.pmf(c, mu_i)
            Prob_Yk += (pi_i * Prob_Yk_i)
        xy.append((c, Prob_Yk))
    (x, y) = zip(*xy)
    asciiplot(x, y, title='Model')