<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def epsilon(data, correction='gg'):
 """Epsilon adjustement factor for repeated measures. Parameters data : pd.DataFrame DataFrame containing the repeated measurements. ``data`` must be in wide-format. To convert from wide to long format, use the :py:func:`pandas.pivot_table` function. correction : string Specify the epsilon version :: 'gg' : Greenhouse-Geisser 'hf' : Huynh-Feldt 'lb' : Lower bound Returns ------- eps : float Epsilon adjustement factor. Notes ----- The **lower bound** for epsilon is: .. math:: lb = \\frac{1}{k - 1} where :math:`k` is the number of groups (= data.shape[1]). The **Greenhouse-Geisser epsilon** is given by: .. math:: \\epsilon_{GG} = \\frac{k^2(\\overline{diag(S)} - \\overline{S})^2} {(k-1)(\\sum_{i=1}^{k}\\sum_{j=1}^{k}s_{ij}^2 - 2k\\sum_{j=1}^{k} \\overline{s_i}^2 + k^2\\overline{S}^2)} where :math:`S` is the covariance matrix, :math:`\\overline{S}` the grandmean of S and :math:`\\overline{diag(S)}` the mean of all the elements on the diagonal of S (i.e. mean of the variances). The **Huynh-Feldt epsilon** is given by: .. math:: \\epsilon_{HF} = \\frac{n(k-1)\\epsilon_{GG}-2}{(k-1) (n-1-(k-1)\\epsilon_{GG})} where :math:`n` is the number of subjects. References .. [1] http://www.real-statistics.com/anova-repeated-measures/sphericity/ Examples -------- 0.5587754577585018 0.6223448311539781 0.5 """

    # Covariance matrix
    S = data.cov()
    n = data.shape[0]
    k = data.shape[1]

    # Lower bound
    if correction == 'lb':
        if S.columns.nlevels == 1:
            return 1 / (k - 1)
        elif S.columns.nlevels == 2:
            ka = S.columns.levels[0].size
            kb = S.columns.levels[1].size
            return 1 / ((ka - 1) * (kb - 1))

    # Compute GGEpsilon
    # - Method 1
    mean_var = np.diag(S).mean()
    S_mean = S.mean().mean()
    ss_mat = (S**2).sum().sum()
    ss_rows = (S.mean(1)**2).sum().sum()
    num = (k * (mean_var - S_mean))**2
    den = (k - 1) * (ss_mat - 2 * k * ss_rows + k**2 * S_mean**2)
    eps = np.min([num / den, 1])
    # - Method 2
    # S_pop = S - S.mean(0)[:, None] - S.mean(1)[None, :] + S.mean()
    # eig = np.linalg.eigvalsh(S_pop)[1:]
    # V = eig.sum()**2 / np.sum(eig**2)
    # eps = V / (k - 1)

    # Huynh-Feldt
    if correction == 'hf':
        num = n * (k - 1) * eps - 2
        den = (k - 1) * (n - 1 - (k - 1) * eps)
        eps = np.min([num / den, 1])
    return eps