def METHOD_NAME(X, y, sample_size, max_iterations, eps=1e-12):
    # randomly draw 2 points for the initial fit
    sample = np.random.choice(len(y), 2, replace=False)
    # perform the initial fit
    beta = linear_regression(X[sample], y[sample])
    if np.isnan(beta[0]):
        return beta, np.nan
    error = residual_sum_of_squares(X[sample], y[sample], beta)
    for i in range(max_iterations):
        squared_residuals = residuals(X, y, beta) ** 2
        # select the subset with the smallest squared residuals
        sample = np.argsort(squared_residuals)[:sample_size]
        # redo regression with new sample
        beta = linear_regression(X[sample], y[sample])
        last_error = error
        error = residual_sum_of_squares(X[sample], y[sample], beta)
        # test for convergences
        if abs(last_error - error) < eps:
            break
    return beta, error