docs: 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