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

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