docs: def normalizedFluctuationCorrelationFunction(A_n, B_n=None, N_max=None, norm=True):
      """Compute the normalized fluctuation (cross) correlation function of (two) stationary timeseries.
  
      C(t) = (<A(t) B(t)> - <A><B>) / (<AB> - <A><B>)
  
      This may be useful in diagnosing odd time-correlations in timeseries data.
  
      Parameters
      ----------
      A_n : np.ndarray
          A_n[n] is nth value of timeseries A.  Length is deduced from vector.
      B_n : np.ndarray
          B_n[n] is nth value of timeseries B.  Length is deduced from vector.
      N_max : int, default=None
          if specified, will only compute correlation function out to time lag of N_max
      norm: bool, optional, default=True
          if False will return the unnormalized correlation function D(t) = <A(t) B(t)>
  
      Returns
      -------
      C_n : np.ndarray
          C_n[n] is the normalized fluctuation auto- or cross-correlation function for timeseries A(t) and B(t).
  
      Notes
      -----
      The same timeseries can be used for both A_n and B_n to get the autocorrelation statistical inefficiency.
      This procedure may be slow.
      The statistical error in C_n[n] will grow with increasing n.  No effort is made here to estimate the uncertainty.
  
      References
      ----------
      [1] J. D. Chodera, W. C. Swope, J. W. Pitera, C. Seok, and K. A. Dill. Use of the weighted
      histogram analysis method for the analysis of simulated and parallel tempering simulations.
      JCTC 3(1):26-41, 2007.
  
      Examples
      --------
  
      Estimate normalized fluctuation correlation function.
  
      >>> from pymbar import testsystems
      >>> A_t = testsystems.correlated_timeseries_example(N=10000, tau=5.0)
      >>> C_t = normalizedFluctuationCorrelationFunction(A_t, N_max=25)
  
      """
  
      # If B_n is not specified, set it to be identical to A_n.
      if B_n is None:
          B_n = A_n
  
      # Create np copies of input arguments.
      A_n = np.array(A_n)
      B_n = np.array(B_n)
  
      # Get the length of the timeseries.
      N = A_n.size
  
      # Set maximum time to compute correlation functon for.
      if (not N_max) or (N_max > N - 1):
          N_max = N - 1
  
      # Be sure A_n and B_n have the same dimensions.
      if(A_n.shape != B_n.shape):
          raise ParameterError('A_n and B_n must have same dimensions.')
  
      # Initialize statistical inefficiency estimate with uncorrelated value.
      g = 1.0
  
      # Compute means and variance.
      mu_A = A_n.mean()
      mu_B = B_n.mean()
  
      # Make temporary copies at high precision with means subtracted off.
      dA_n = A_n.astype(np.float64) - mu_A
      dB_n = B_n.astype(np.float64) - mu_B
  
      # sigma2_AB = sum((A_n-mu_A) * (B_n-mu_B)) / (float(N)-1.0) # unbiased estimator
      sigma2_AB = (dA_n * dB_n).mean()  # standard estimator to ensure C(0) = 1
      if(sigma2_AB == 0):
          raise ParameterError('Sample covariance sigma_AB^2 = 0 -- cannot compute statistical inefficiency')
  
      # allocate storage for normalized fluctuation correlation function
      C_n = np.zeros([N_max + 1], np.float64)
  
      # Compute normalized correlation function.
      t = 0
      for t in range(0, N_max + 1):
          # compute normalized fluctuation correlation function at time t
          C_n[t] = np.sum(dA_n[0:(N - t)] * dB_n[t:N] + dB_n[0:(N - t)] * dA_n[t:N]) / (2.0 * float(N - t) * sigma2_AB)
  
      # Return the computed correlation function
      if norm:
          return C_n
      else:
          return C_n*sigma2_AB + mu_A*mu_B