docs: def direct_age_standardization(e, b, s, n, alpha=0.05):
      """A utility function to compute rate through direct age standardization
  
      Parameters
      ----------
      e          : array
                   (n*h, 1), event variable measured for each age group across n spatial units
      b          : array
                   (n*h, 1), population at risk variable measured for each age group across n spatial units
      s          : array
                   (n*h, 1), standard population for each age group across n spatial units
      n          : integer
                   the number of spatial units
      alpha      : float
                   significance level for confidence interval
  
      Notes
      -----
      e, b, and s are arranged in the same order
  
      Returns
      -------
      list
          a list of n tuples; a tuple has a rate and its lower and upper limits
          age standardized rates and confidence intervals
  
      Examples
      --------
  
      Creating an array of an event variable (e.g., the number of cancer patients)
      for 2 regions in each of which 4 age groups are available.
      The first 4 values are event values for 4 age groups in the region 1,
      and the next 4 values are for 4 age groups in the region 2.
  
      >>> e = np.array([30, 25, 25, 15, 33, 21, 30, 20])
  
      Creating another array of a population-at-risk variable (e.g., total population)
      for the same two regions.
      The order for entering values is the same as the case of e.
  
      >>> b = np.array([1000, 1000, 1100, 900, 1000, 900, 1100, 900])
  
      For direct age standardization, we also need the data for standard population.
      Standard population is a reference population-at-risk (e.g., population distribution for the U.S.)
      whose age distribution can be used as a benchmarking point for comparing age distributions
      across regions (e.g., population distribution for Arizona and California).
      Another array including standard population is created.
  
      >>> s = np.array([1000, 900, 1000, 900, 1000, 900, 1000, 900])
  
      Specifying the number of regions.
  
      >>> n = 2
  
      Applying direct_age_standardization function to e and b
  
      >>> a, b = [i[0] for i in direct_age_standardization(e, b, s, n)]
      >>> round(a, 4)
      0.0237
      >>> round(b, 4)
      0.0267
  
      """
      age_weight = (1.0 / b) * (s * 1.0 / sum_by_n(s, 1.0, n).repeat(len(s) // n))
      adjusted_r = sum_by_n(e, age_weight, n)
      var_estimate = sum_by_n(e, np.square(age_weight), n)
      g_a = np.square(adjusted_r) / var_estimate
      g_b = var_estimate / adjusted_r
      k = [age_weight[i:i + len(b) // n].max() for i in range(0, len(b),
                                                             len(b) // n)]
      g_a_k = np.square(adjusted_r + k) / (var_estimate + np.square(k))
      g_b_k = (var_estimate + np.square(k)) / (adjusted_r + k)
      summed_b = sum_by_n(b, 1.0, n)
      res = []
      for i in range(len(adjusted_r)):
          if adjusted_r[i] == 0:
              upper = 0.5 * chi2.ppf(1 - 0.5 * alpha)
              lower = 0.0
          else:
              lower = gamma.ppf(0.5 * alpha, g_a[i], scale=g_b[i])
              upper = gamma.ppf(1 - 0.5 * alpha, g_a_k[i], scale=g_b_k[i])
          res.append((adjusted_r[i], lower, upper))
      return res