docs: def draw_noisy_time_series(self, SNR=1.0, red_noise_ratio=0.25, outlier_ratio=0.0):
          """
          A function to draw a noisy time series based on the clean model 
          such that 
              y_noisy = y + yw + yr,
          where yw is white noise, yr is red noise and y will be rescaled 
          so that y_noisy complies with the specified signal-to-noise 
          ratio (SNR). 
          
          Parameters
          ---------
          SNR: float
              Signal-to-noise ratio of the resulting contaminated signal in
              decibels [dB]. SNR is defined as 
                  SNR = 10*log(var_signal/var_noise), hence
              NR var_signal/var_noise
                   10  10
                    7   5
                    3   2
                    0   1
                   -3   0.5
                   -7   0.2
                  -10   0.1
          red_noise_variance: float in [0, 1]
              The variance of the red noise component is set according to 
              Var(yw)*red_noise_ratio. Set this to zero to obtain uncertainties
              that explain the noise perfectly
          outlier_ratio: float in [0, 1]
              Percentage of outlier data points
              
          Returns
          -------
          t: ndarray
              Vector containing the time instants
          y_noisy: ndarray
              Vector containing the contaminated signal
          s: ndarray
              Vector containing the uncertainties associated to the white 
              noise component
  
          """
          if outlier_ratio < 0.0 or outlier_ratio > 1.0:
              raise ValueError("Outlier ratio must be in [0, 1]")
          if red_noise_ratio < 0.0:
              raise ValueError("Red noise ratio must be positive")
          np.random.seed(self.rseed)
          t = self.t
          y_clean = self.y_clean
          N = len(t)
          # First we generate s 
          s, mean_s_squared = generate_uncertainties(N, rseed=self.rseed)
          #print(mean_s_squared)
          #print(np.mean(s**2))
          # Draw a heteroscedastic white noise vector
          white_noise = np.random.multivariate_normal(np.zeros(N,), np.diag(s**2))
          # Now we generate a colored noise vector which is unaccounted by s
          red_noise_variance = mean_s_squared*red_noise_ratio
          # First order markovian process to generate 
          red_noise = first_order_markov_process(t, red_noise_variance, 1.0, rseed=self.rseed)
          
          # The following is not ok for irregularly sampled time series because
          # it assumes constant dt=1
          #phi=0.5
          #red_noise = np.random.randn(N)*np.sqrt(red_noise_variance)
          #for i in range(1, N):
          #    red_noise[i] = phi*red_noise[i-1] + np.sqrt(1 - phi**2)*red_noise[i]
          
          # The final noise vector
          #print("%f %f" % (np.var(white_noise)*red_noise_ratio, np.var(red_noise)))
          noise = white_noise + red_noise 
          var_noise = mean_s_squared + red_noise_variance
          SNR_unitless = 10.0**(SNR/10.0)
          self.A = np.sqrt(SNR_unitless*var_noise)
          y = self.A*y_clean
          y_noisy = y + noise
          # Add outliers with a certain percentage
          rperm = np.where(np.random.uniform(size=N) < outlier_ratio)[0]    
          outlier = np.random.uniform(5.0*np.std(y), 10.0*np.std(y), size=len(rperm))
          y_noisy[rperm] += outlier
          return t, y_noisy, s