Divide two variables with asymmetric uncertainties

rel_asym_err.py

import scipy.optimize

def mlnL_Barlow(x, x_mean, sigma_p, sigma_n):
    # from arXiv:physics/0406120
    return 0.5 * ((x - x_mean)**2) / (sigma_p * sigma_n + (sigma_p - sigma_n) * (x - x_mean))

assert mlnL_Barlow(7+5, 7, 5, 3) == 0.5
assert mlnL_Barlow(7-3, 7, 5, 3) == 0.5


def find_rel(x_mean, x_sigma_p, x_sigma_n, y_mean, y_sigma_p, y_sigma_n):
    def mlnL_rel(r):
        def f(y):
            return mlnL_Barlow(r*y, x_mean, x_sigma_p, x_sigma_n) + mlnL_Barlow(y, y_mean, y_sigma_p, y_sigma_n)
        res = scipy.optimize.minimize(f, x_mean/r, method='nelder-mead')
        return res.fun
    sol_p = scipy.optimize.newton_krylov(lambda x: (mlnL_rel(x) - 0.5), (x_mean + x_sigma_p) / y_mean)
    sol_n = scipy.optimize.newton_krylov(lambda x: (mlnL_rel(x) - 0.5), (x_mean - x_sigma_n) / y_mean)
    assert sol_p > sol_n

    res = scipy.optimize.minimize(mlnL_rel, x_mean / y_mean, method='nelder-mead')
    r = x_mean / y_mean
    assert abs(res.x - r) < 0.0001
    assert sol_p > r
    assert sol_n < r
    return (r, sol_p - r, r - sol_n)

    from matplotlib import pyplot
    import numpy as np
    rs = np.linspace(sol_n, sol_p, 1000)
    pyplot.plot(rs, map(lambda r: mlnL_rel(r), rs))
    pyplot.show()