def METHOD_NAME(y, X=None, max_order=6):
    """Polynomial Regression.
    Find the optimal order for polynomial fitting. Currently, the only method implemented is
    RMSE minimization.
    Parameters
    ----------
    y : Union[list, np.array, pd.Series]
        The response variable (the y axis).
    X : Union[list, np.array, pd.Series]
        Explanatory variable (the x axis). If 'None', will treat y as a continuous signal.
    max_order : int
        The maximum order to test.
    Returns
    -------
    int
        Optimal order.
    See Also
    ----------
    fit_polynomial
    Examples
    ---------
      import neurokit2 as nk
      y = np.cos(np.linspace(start=0, stop=10, num=100))
      nk.fit_polynomial_findorder(y, max_order=10)
    9
    """
    # TODO: add cross-validation or some kind of penalty to prevent over-fitting?
    if X is None:
        X = np.linspace(0, 100, len(y))
    best_rmse = 0
    for order in range(max_order):
        y_predicted, _ = _fit_polynomial(y, X, order=order)
        rmse = fit_rmse(y, y_predicted)
        if rmse < best_rmse or best_rmse == 0:
            best_order = order
    return best_order