code: def rs(data, n, unbiased=True):
    
    data = np.asarray(data)
    total_N = len(data)
    m = total_N // n # number of sequences
    # cut values at the end of data to make the array divisible by n
    data = data[:total_N - (total_N % n)]
    # split remaining data into subsequences of length n
    seqs = np.reshape(data, (m, n))
    # calculate means of subsequences
    means = np.mean(seqs, axis=1)
    # normalize subsequences by substracting mean
    y = seqs - means.reshape((m, 1))
    # build cumulative sum of subsequences
    y = np.cumsum(y, axis=1)
    # find ranges
    r = np.max(y, axis=1) - np.min(y, axis=1)
    # find standard deviation
    # we should use the unbiased estimator, since we do not know the true mean
    s = np.std(seqs, axis=1, ddof=1 if unbiased else 0)
    # some ranges may be zero and have to be excluded from the analysis
    idx = np.where(r != 0)
    r = r[idx]
    s = s[idx]
    # it may happen that all ranges are zero (if all values in data are equal)
    if len(r) == 0:
      return np.nan
    else:
      # return mean of r/s along subsequence index
      return np.mean(r / s)
docs: Calculates an individual R/S value in the rescaled range approach for
    a given n.
  
    Note: This is just a helper function for hurst_rs and should not be called
    directly.
  
    Args:
      data (array-like of float):
        time series
      n (float):
        size of the subseries in which data should be split
  
    Kwargs:
      unbiased (boolean):
        if True, the standard deviation based on the unbiased variance
        (1/(N-1) instead of 1/N) will be used. This should be the default choice,
        since the true mean of the sequences is not known. This parameter should
        only be changed to recreate results of other implementations.
  
    Returns:
      float:
        (R/S)_n