Error in constraint which makes lower bound higher than upper bound

[Ignatiocalvin]

Hello,

I have encountered an error within the constraint
max(row['lb-'+sam] * 2, 0), min(1, 2 * row['ub-'+sam])])

In line 13 of process_A7_new.ipynb within the binomial_hdpr() which prevents ./generatemutationtrees from being run.

I have included the dataset to replicate this error as well as some minor changes to the function to be able to run them without errors. (encountered errors were floating point error, error when ref- & alt-counts were 0 and string index out of range errors).

Specifically:

def binomial_hpdr()
...
if N != 0:
        mode = (n+a-1.)/(N+a+b-2.)
    else: return("")
...
 p_range = numpy.linspace(min_p, max_p, int(n_pbins+1))

and the conditions for the get_ functions

Here are the changes that have been done:

def binomial_hpdr(n, N, pct, a=1, b=1, n_pbins=1e3):
    """
    Function computes the posterior mode along with the upper and lower bounds of the
    **Highest Posterior Density Region**.

    Parameters
    ----------
    n: number of successes 
    N: sample size 
    pct: the size of the confidence interval (between 0 and 1)
    a: the alpha hyper-parameter for the Beta distribution used as a prior (Default=1)
    b: the beta hyper-parameter for the Beta distribution used as a prior (Default=1)
    n_pbins: the number of bins to segment the p_range into (Default=1e3)

    Returns
    -------
    A tuple that contains the mode as well as the lower and upper bounds of the interval
    (mode, lower, upper)

    """
    # fixed random variable object for posterior Beta distribution
    rv = beta(n+a, N-n+b)
    # determine the mode and standard deviation of the posterior
    stdev = rv.stats('v')**0.5
    if N != 0:
        mode = (n+a-1.)/(N+a+b-2.)
    else: return("")
    # compute the number of sigma that corresponds to this confidence
    # this is used to set the rough range of possible success probabilities
    n_sigma = numpy.ceil(norm.ppf( (1+pct)/2. ))+1
    # set the min and max values for success probability 
    max_p = mode + n_sigma * stdev
    if max_p > 1:
        max_p = 1.
    min_p = mode - n_sigma * stdev
    if min_p > 1:
        min_p = 1.
    # make the range of success probabilities
    p_range = numpy.linspace(min_p, max_p, int(n_pbins+1))
    # construct the probability mass function over the given range
    if mode > 0.5:
        sf = rv.sf(p_range)
        pmf = sf[:-1] - sf[1:]
    else:
        cdf = rv.cdf(p_range)
        pmf = cdf[1:] - cdf[:-1]
    # find the upper and lower bounds of the interval 
    sorted_idxs = numpy.argsort( pmf )[::-1]
    cumsum = numpy.cumsum( numpy.sort(pmf)[::-1] )
    j = numpy.argmin( numpy.abs(cumsum - pct) )
    upper = p_range[ (sorted_idxs[:j+1]).max()+1 ]
    lower = p_range[ (sorted_idxs[:j+1]).min() ]    

    return (mode, lower, upper)

def get_ub(row):
    v=binomial_hpdr(row["alt_counts"], row["alt_counts"]+row["ref_counts"], corrected_confidence)
    if v =='':
        return 0
    else:
        return v[2]
    

def get_lb(row):
    v=binomial_hpdr(row["alt_counts"], row["alt_counts"]+row["ref_counts"], corrected_confidence)
    if v =='':
        return 0
    else:
        mval = v[1]
        if mval < 0.01: mval = 0
        return mval

def get_mean(row):
    v=binomial_hpdr(row["alt_counts"], row["alt_counts"]+row["ref_counts"], corrected_confidence)
    if v =='':
        return 0
    else:
        mval = v[0]
        return mval

Your help and insight would be very much appreciated and I would like to thank you for your time.

test.csv