Note that this function requires Scipy.

def sample_double_norm(mean, std_upper, std_lower, size):
    """Note that this function requires Scipy."""
    from scipy.special import erfinv

    # There's probably a better way to do this. We first draw percentiles
    # uniformly between 0 and 1. We want the peak of the distribution to occur
    # at `mean`. However, if we assign 50% of the samples to the lower half
    # and 50% to the upper half, the side with the smaller variance will be
    # overrepresented because of the 1/sigma normalization of the Gaussian
    # PDF. Therefore we need to divide points between the two halves with a
    # fraction `cutoff` (defined below) going to the lower half. Having
    # partitioned them this way, we can then use the standard Gaussian
    # quantile function to go from percentiles to sample values -- except that
    # we must remap from [0, cutoff] to [0, 0.5] and from [cutoff, 1] to [0.5,
    # 1].

    samples = np.empty(size)
    percentiles = np.random.uniform(0., 1., size)
    cutoff = std_lower / (std_lower + std_upper)

    w = (percentiles < cutoff)
    percentiles[w] *= 0.5 / cutoff
    samples[w] = mean + np.sqrt(2) * std_lower * erfinv(2 * percentiles[w] - 1)

    w = ~w
    percentiles[w] = 1 - (1 - percentiles[w]) * 0.5 / (1 - cutoff)
    samples[w] = mean + np.sqrt(2) * std_upper * erfinv(2 * percentiles[w] - 1)

    return samples