docs: def power_law_noise(t, variance):
      """
      CAUTION: NOT THOROUGHLY TESTED
      Generates a red noise vector by first generating a covariance 
      function based on the expected red noise spectra. Red noise spectrum
      follows a power law
      
          S(f) = A*f^(-2),
      
      where A is a constant. Pink noise can be obtained if the exponent is
      changed to -1, and white noise if it is changed to 0.
      
      This Power Spectral density (PSD) is real and even then
      
          r(tau) = AT sum cos(2pi k tau/T)/k**2,
      
      where we assume f[k] = k*df, df = 1/T, Fs = N/T.
      
      The basel problem gives us: sum (1/k**2) approx pi**2/6, hence we can
      set c according to the desired variance at lag=0
       
      After the covariance is computed we can draw from a multivariate
      normal distribution to obtain a noise vector
      
      Parameters
      ---------
      t: ndarray
          A time vector for which the red noise vector will be sampled
      variance: positive float
          variance of the resulting red noise vector
          
      Returns
      -------
      red_noise: ndarray
          Vector containing the red noise realizations
          
      See also
      --------
      first_order_markov_process
      
      """
      if variance < 0.0:
          raise ValueError("Variance must be positive")
      mu = np.zeros(shape=(N,))
      if variance == 0.0:
          return mu
      N = len(t)    
      T = (t[-1] - t[0])
      c = (6.0*variance)/(T*np.pi**2)
      f = np.arange(1.0/T, 0.5*N/T, step=1.0/T)  # We ommit f=0.0
      k = f*T
      dt = np.repeat(np.reshape(t, (1, -1)), N, axis=0)
      dt = np.absolute(dt - dt.T)  # This is NxN
      S = np.zeros(shape=(N, N))
      for i in range(0, N):
          for j in range(i, N):
              S[i, j] = np.sum(np.cos(2.0*np.pi*k*dt[i, j]/T)/k**2)*T*c
              S[j, i] = S[i, j] 
      red_noise = np.random.multivariate_normal(mu, S)
      return red_noise