docs: def indirect_age_standardization(e, b, s_e, s_b, n, alpha=0.05):
      """A utility function to compute rate through indirect 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_e        : array
                   (n*h, 1), event variable measured for each age group across n spatial units in a standard population
      s_b        : array
                   (n*h, 1), population variable measured for each age group across n spatial units in a standard population
      n          : integer
                   the number of spatial units
      alpha      : float
                   significance level for confidence interval
  
      Notes
      -----
      e, b, s_e, and s_b 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 rate
  
      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([100, 100, 110, 90, 100, 90, 110, 90])
  
      For indirect age standardization, we also need the data for standard population and event.
      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., popoulation distribution for Arizona and California).
      When the same concept is applied to the event variable,
      we call it standard event (e.g., the number of cancer patients in the U.S.).
      Two additional arrays including standard population and event are created.
  
      >>> s_e = np.array([100, 45, 120, 100, 50, 30, 200, 80])
      >>> s_b = np.array([1000, 900, 1000, 900, 1000, 900, 1000, 900])
  
      Specifying the number of regions.
  
      >>> n = 2
  
      Applying indirect_age_standardization function to e and b
  
      >>> [i[0] for i in indirect_age_standardization(e, b, s_e, s_b, n)]
      [0.23723821989528798, 0.2610803324099723]
  
      """
      smr = standardized_mortality_ratio(e, b, s_e, s_b, n)
      s_r_all = sum(s_e * 1.0) / sum(s_b * 1.0)
      adjusted_r = s_r_all * smr
  
      e_by_n = sum_by_n(e, 1.0, n)
      log_smr = np.log(smr)
      log_smr_sd = 1.0 / np.sqrt(e_by_n)
      norm_thres = norm.ppf(1 - 0.5 * alpha)
      log_smr_lower = log_smr - norm_thres * log_smr_sd
      log_smr_upper = log_smr + norm_thres * log_smr_sd
      smr_lower = np.exp(log_smr_lower) * s_r_all
      smr_upper = np.exp(log_smr_upper) * s_r_all
      res = list(zip(adjusted_r, smr_lower, smr_upper))
      return res