<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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. 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. 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. Specifying the number of regions. Applying direct_age_standardization function to e and b 0.0237 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