Как найти наилучшую степень многочлена?

Вы должны предоставить данные для X/Y в следующий раз или что-то фиктивное, это будет быстрее и предоставит вам конкретное решение. На данный момент я создал фиктивное уравнение формы y = X**4 + X**3 + X + 1.

Есть много способов улучшить это, но быстрая итерация для поиска наилучшей степени заключается в том, чтобы просто подобрать ваши данные для каждой степени и выбрать степень с наилучшей производительностью (например, с самым низким RMSE).

Вы также можете поиграть с тем, как вы решите хранить данные обучения/тестирования/проверки.

import numpy as np
import matplotlib.pyplot as plt 

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

X = np.arange(100).reshape(100, 1)
y = X**4 + X**3 + X + 1

x_train, x_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

rmses = []
degrees = np.arange(1, 10)
min_rmse, min_deg = 1e10, 0

for deg in degrees:

    # Train features
    poly_features = PolynomialFeatures(degree=deg, include_bias=False)
    x_poly_train = poly_features.fit_transform(x_train)

    # Linear regression
    poly_reg = LinearRegression()
    poly_reg.fit(x_poly_train, y_train)

    # Compare with test data
    x_poly_test = poly_features.fit_transform(x_test)
    poly_predict = poly_reg.predict(x_poly_test)
    poly_mse = mean_squared_error(y_test, poly_predict)
    poly_rmse = np.sqrt(poly_mse)
    rmses.append(poly_rmse)
    
    # Cross-validation of degree
    if min_rmse > poly_rmse:
        min_rmse = poly_rmse
        min_deg = deg

# Plot and present results
print('Best degree {} with RMSE {}'.format(min_deg, min_rmse))
        
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(degrees, rmses)
ax.set_yscale('log')
ax.set_xlabel('Degree')
ax.set_ylabel('RMSE')

Это напечатает:

Лучшая степень 4 с RMSE 1.27689038706e-08

В качестве альтернативы вы также можете создать новый класс, который выполняет подбор полиномов, и передать его в GridSearchCV с набором параметров.

На мой взгляд, лучший способ найти оптимальную степень подгонки кривой или вообще подходящую модель — использовать модуль GridSearchCV из библиотеки scikit-learn.

Вот пример использования этой библиотеки:

Во-первых, давайте определим метод для выборки случайных данных:

def make_data(N, err=1.0, rseed=1):

    rng = np.random.RandomState(rseed)
    X = rng.rand(N, 1) ** 2
    y = 1. / (X.ravel() + 0.3)
    if err > 0:
        y += err * rng.randn(N)
    return X, y

Построить трубопровод:

def PolynomialRegression(degree=2, **kwargs):
    return make_pipeline(PolynomialFeatures(degree), LinearRegression(**kwargs))

Создайте данные и вектор (X_test) для тестирования и визуализации:

X, y = make_data(200)
X_test = np.linspace(-0.1, 1.1, 200)[:, None]

Определите параметры GridSearchCV:

param_grid = {'polynomialfeatures__degree': np.arange(20),
'linearregression__fit_intercept': [True, False],
'linearregression__normalize': [True, False]}
grid = GridSearchCV(PolynomialRegression(), param_grid, cv=7)
grid.fit(X, y)

Получите лучшие параметры от нашей модели:

model = grid.best_estimator_
model

Pipeline(memory=None,
     steps=[('polynomialfeatures', PolynomialFeatures(degree=4, include_bias=True, interaction_only=False)), ('linearregression', LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False))])

Подберите модель с данными X и y и используйте вектор, чтобы предсказать значения:

y_test = model.fit(X, y).predict(X_test)

Визуализируйте результат:

plt.scatter(X, y)
plt.plot(X_test.ravel(), y_test, 'r')

Наилучший результат

Полный фрагмент кода:

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import GridSearchCV

def make_data(N, err=1.0, rseed=1):

    rng = np.random.RandomState(rseed)
    X = rng.rand(N, 1) ** 2
    y = 1. / (X.ravel() + 0.3)
    if err > 0:
        y += err * rng.randn(N)
    return X, y

def PolynomialRegression(degree=2, **kwargs):
    return make_pipeline(PolynomialFeatures(degree), LinearRegression(**kwargs))


X, y = make_data(200)
X_test = np.linspace(-0.1, 1.1, 200)[:, None]

param_grid = {'polynomialfeatures__degree': np.arange(20),
'linearregression__fit_intercept': [True, False],
'linearregression__normalize': [True, False]}
grid = GridSearchCV(PolynomialRegression(), param_grid, cv=7)
grid.fit(X, y)

model = grid.best_estimator_

y_test = model.fit(X, y).predict(X_test)

plt.scatter(X, y)
plt.plot(X_test.ravel(), y_test, 'r')

Вот здесь-то и вступает в дело выбор байесовской модели. Это дает вам наиболее вероятную модель, учитывая как сложность модели, так и соответствие данных. Я очень устал, поэтому быстрый ответ - использовать BIC (байесовский информационный критерий):

k = number of variables in the model
n = number of observations
sse = sum(residuals**2)
BIC = n*ln(sse/n) + k*ln(n) 

Этот BIC (или AIC и т. д.) даст вам лучшую модель