docs: def fit_creatine(self, reject_outliers=3.0, fit_lb=2.7, fit_ub=3.5):
          """
          Fit a model to the portion of the summed spectra containing the
          creatine and choline signals.
  
          Parameters
          ----------
          reject_outliers : float or bool
             If set to a float, this is the z score threshold for rejection (on
             any of the parameters). If set to False, no outlier rejection
  
          fit_lb, fit_ub : float
             What part of the spectrum (in ppm) contains the creatine peak.
             Default (2.7, 3.5)
  
          Note
          ----
          We use upper and lower bounds that are a variation on the bounds
          mentioned on the GANNET ISMRM2013 poster [1]_.
  
          [1] RAE Edden et al (2013). Gannet GABA analysis toolkit. ISMRM
          conference poster.
  
          """
          # We fit a two-lorentz function to this entire chunk of the spectrum,
          # to catch both choline and creatine
          model, signal, params = ana.fit_two_lorentzian(self.sum_spectra,
                                                         self.f_ppm,
                                                         lb=fit_lb,
                                                         ub=fit_ub)
  
          # Use an array of ones to index everything but the outliers and nans:
          ii = np.ones(signal.shape[0], dtype=bool)
          # Reject outliers:
          if reject_outliers:
              model, signal, params, ii = self._outlier_rejection(params,
                                                                  model,
                                                                  signal,
                                                                  ii)
              
          # We'll keep around a private attribute to tell us which transients
          # were good (this is for both creatine and choline):
          self._cr_transients = np.where(ii)
          
          # Now we separate choline and creatine params from each other (remember
          # that they both share offset and drift!):
          self.choline_params = params[:, (0,2,4,6,8,9)]
          self.creatine_params = params[:, (1,3,5,7,8,9)]
          
          self.cr_idx = ut.make_idx(self.f_ppm, fit_lb, fit_ub)
  
          # We'll need to generate the model predictions from these parameters,
          # because what we're holding in 'model' is for both together:
          self.choline_model = np.zeros((self.creatine_params.shape[0],
                                      np.abs(self.cr_idx.stop-self.cr_idx.start)))
  
          self.creatine_model = np.zeros((self.choline_params.shape[0],
                                      np.abs(self.cr_idx.stop-self.cr_idx.start)))
          
          for idx in range(self.creatine_params.shape[0]):
              self.creatine_model[idx] = ut.lorentzian(self.f_ppm[self.cr_idx],
                                                       *self.creatine_params[idx])
              self.choline_model[idx] = ut.lorentzian(self.f_ppm[self.cr_idx],
                                                      *self.choline_params[idx])
          self.creatine_signal = signal
          self.creatine_auc = self._calc_auc(ut.lorentzian,
                                             self.creatine_params,
                                             self.cr_idx)
          self.choline_auc = self._calc_auc(ut.lorentzian,
                                            self.choline_params,
                                            self.cr_idx)