docs: def normalizedFluctuationCorrelationFunctionMultiple(A_kn, B_kn=None, N_max=None, norm=True, truncate=False):
      """Compute the normalized fluctuation (cross) correlation function of (two) timeseries from multiple timeseries samples.
  
      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_kn : Python list of numpy arrays
          A_kn[k] is the kth timeseries, and A_kn[k][n] is nth value of timeseries k.  Length is deduced from arrays.
      B_kn : Python list of numpy arrays
          B_kn[k] is the kth timeseries, and B_kn[k][n] is nth value of timeseries k.  B_kn[k] must have same length as A_kn[k]
      N_max : int, optional, 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 unnormalized D(t) = <A(t) B(t)>
      truncate: bool, optional, default=False
          if True, will stop calculating the correlation function when it goes below 0
  
      Returns
      -------
      C_n[n] : np.ndarray
          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 a portion of the normalized fluctuation autocorrelation function from multiple timeseries of different length.
  
      >>> from pymbar import testsystems
      >>> N_k = [1000, 2000, 3000, 4000, 5000]
      >>> tau = 5.0 # exponential relaxation time
      >>> A_kn = [ testsystems.correlated_timeseries_example(N=N, tau=tau) for N in N_k ]
      >>> C_n = normalizedFluctuationCorrelationFunctionMultiple(A_kn, N_max=25)
  
      """
  
      # If B_kn is not specified, define it to be identical with A_kn.
      if B_kn is None:
          B_kn = A_kn
  
      # TODO: Change this to support other iterable types, like sets.
      # Make sure A_kn and B_kn are both lists
      if (type(A_kn) is not list) or (type(B_kn) is not list):
          raise ParameterError("A_kn and B_kn must each be a list of numpy arrays.")
  
      # Ensure the same number of timeseries are stored in A_kn and B_kn.
      if (len(A_kn) != len(B_kn)):
          raise ParameterError("A_kn and B_kn must contain corresponding timeseries -- different numbers of timeseries detected in each.")
  
      # Determine number of timeseries stored.
      K = len(A_kn)
  
      # Ensure both observable trajectories in each timeseries are of the same length.
      for k in range(K):
          A_n = A_kn[k]
          B_n = B_kn[k]
          if A_n.size != B_n.size:
              raise ParameterError("A_kn and B_kn must contain corresponding timeseries -- lack of correspondence in timeseries lenghts detected.")
  
      # Get the length of each timeseries.
      N_k = np.zeros([K], np.int32)
      for k in range(K):
          N_k[k] = A_kn[k].size
  
      # Determine total number of samples.
      N = np.sum(N_k)
  
      # Set maximum time to compute correlation functon for.
      if (not N_max) or (N_max > max(N_k) - 1):
          N_max = max(N_k) - 1
  
      # Compute means.
      mu_A = 0.0
      mu_B = 0.0
      for k in range(K):
          mu_A += np.sum(A_kn[k])
          mu_B += np.sum(B_kn[k])
      mu_A /= float(N)
      mu_B /= float(N)
  
      # Compute fluctuation timeseries.
      dA_kn = list()
      dB_kn = list()
      for k in range(K):
          dA_n = A_kn[k] - mu_A
          dB_n = B_kn[k] - mu_B
          dA_kn.append(dA_n)
          dB_kn.append(dB_n)
  
      # Compute covariance.
      sigma2_AB = 0.0
      for k in range(K):
          sigma2_AB += np.sum(dA_kn[k] * dB_kn[k])
      sigma2_AB /= float(N)
  
      # allocate storage for normalized fluctuation correlation function
      C_n = np.zeros([N_max + 1], np.float64)
  
      # Accumulate the integrated correlation time by computing the normalized correlation time at
      # increasing values of t.  Stop accumulating if the correlation function goes negative, since
      # this is unlikely to occur unless the correlation function has decayed to the point where it
      # is dominated by noise and indistinguishable from zero.
      t = 0
      negative = False
      for t in range(0, N_max + 1):
          # compute unnormalized correlation function
          numerator = 0.0
          denominator = 0.0
          for k in range(K):
              if (t >= N_k[k]):
                  continue  # skip this trajectory if t is longer than the timeseries
              numerator += np.sum(dA_kn[k][0:(N_k[k] - t)] * dB_kn[k][t:N_k[k]])
              denominator += float(N_k[k] - t)
              if truncate and numerator < 0:
                  negative = True
          C = numerator / denominator
  
          # compute normalized fluctuation correlation function at time t
          C /= sigma2_AB
  
          # Store correlation function.
          C_n[t] = C
  
          if negative:
              break
  
      # Return the computed fluctuation correlation function.
      if norm:
          return C_n[:t]
      else:
          return C_n[:t]*sigma2_AB + mu_A*mu_B