docs: def _gamma_difference_hrf(tr, oversampling=50, time_length=32., onset=0.,
                            delay=6, undershoot=16., dispersion=1.,
                            u_dispersion=1., ratio=0.167):
      """ Compute an hrf as the difference of two gamma functions
  
      Parameters
      ----------
  
      tr : float
          scan repeat time, in seconds
  
      oversampling : int, optional (default=16)
          temporal oversampling factor
  
      time_length : float, optional (default=32)
          hrf kernel length, in seconds
  
      onset: float
          onset time of the hrf
  
      delay: float, optional
          delay parameter of the hrf (in s.)
  
      undershoot: float, optional
          undershoot parameter of the hrf (in s.)
  
      dispersion : float, optional
          dispersion parameter for the first gamma function
  
      u_dispersion : float, optional
          dispersion parameter for the second gamma function
  
      ratio : float, optional
          ratio of the two gamma components
  
      Returns
      -------
      hrf : array of shape(length / tr * oversampling, dtype=float)
           hrf sampling on the oversampled time grid
      """
      from scipy.stats import gamma
      dt = tr / oversampling
      time_stamps = np.linspace(0, time_length, np.rint(float(time_length) / dt).astype(np.int))
      time_stamps -= onset
      hrf = gamma.pdf(time_stamps, delay / dispersion, dt / dispersion) -\
          ratio * gamma.pdf(
          time_stamps, undershoot / u_dispersion, dt / u_dispersion)
      hrf /= hrf.sum()
      return hrf