Given a modal value and a standard deviation, compute corresponding
    parameters for the gamma distribution.

    Intended to be used to replace normal distributions when the value must be
    positive and the uncertainty is comparable to the best value. Conversion
    equations determined from the relations given in the sample_gamma()
    docs.

def find_gamma_params(mode, std):
    """Given a modal value and a standard deviation, compute corresponding
    parameters for the gamma distribution.

    Intended to be used to replace normal distributions when the value must be
    positive and the uncertainty is comparable to the best value. Conversion
    equations determined from the relations given in the sample_gamma()
    docs.

    """
    if mode < 0:
        raise ValueError('input mode must be positive for gamma; got %e' % mode)

    var = std**2
    beta = (mode + np.sqrt(mode**2 + 4 * var)) / (2 * var)
    j = 2 * var / mode**2
    alpha = (j + 1 + np.sqrt(2 * j + 1)) / j

    if alpha <= 1:
        raise ValueError('couldn\'t compute self-consistent gamma parameters: '
                         'mode=%e std=%e alpha=%e beta=%e' % (mode, std, alpha, beta))

    return alpha, beta