Практика 13 · Зміщення проти дисперсії¶

📖 Теорія: відкрий lecture.html у цій же теці. 📝 Домашнє завдання: homework.md. 🧠 Тест: quiz.html.

Лекція каже: помилка розкладається на три доданки — шум, зміщення² і дисперсію, і цей розклад точний. Тут ми його порахуємо руками й переконаємось, що доданки справді сумуються в загальну помилку.

Що зробимо:

  1. Створимо «сорок паралельних світів» — багато навчальних вибірок з одного джерела
  2. Побудуємо віяло прогнозів і побачимо дисперсію очима
  3. Порахуємо bias², variance і шум числом та перевіримо, що вони складаються в очікувану помилку
  4. Побудуємо криві train/test від складності — те, що бачить практик з однією вибіркою
  5. Зробимо те саме через бутстреп, коли вибірка одна-єдина
  6. Покажемо, що дані б'ють лише по дисперсії, і то як $1/n$
  7. Обміняємо дисперсію на зміщення регуляризацією — і виграємо в сумі
  8. Побудуємо те саме готовими інструментами learning_curve і validation_curve

0. Істина, шум і вибірка¶

Влаштовуємо світ, у якому ми знаємо істину. У житті так не буває — саме тому bias і variance неможливо поміряти на реальних даних. Тут можна, бо істину задали ми самі.

$$y = f(x) + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0,\ \sigma^2)$$

Функція $f$ навмисно не поліном: жоден скінченний степінь не відтворить її точно, тому зміщення буде справжнім, а не штучним.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import legendre

NOISE_SIGMA = 0.3          # стандартне відхилення шуму
SAMPLE_SIZE = 40           # скільки точок у кожній навчальній вибірці


def true_function(x):
    """Істина, якої модель не знає. Синусоїда — не поліном, і в цьому вся сіль."""
    return np.sin(1.7 * np.pi * x) + 0.6 * x ** 2


# сітка, на якій міряємо помилку; краї трохи обрізані — там будь-який поліном шаліє
x_test = np.linspace(-0.85, 0.85, 171)
f_test = true_function(x_test)


def draw_training_sample(rng, size=SAMPLE_SIZE):
    """Одна навчальна вибірка: випадкові точки плюс шум."""
    x = rng.uniform(-1, 1, size)
    y = true_function(x) + rng.normal(0, NOISE_SIGMA, size)
    return x, y


demo_rng = np.random.default_rng(1)
x_demo, y_demo = draw_training_sample(demo_rng)

print(f"розмір вибірки: {SAMPLE_SIZE} точок")
print(f"шум: σ = {NOISE_SIGMA}, отже незвідна частина помилки σ² = {NOISE_SIGMA ** 2:.4f}")
print(f"нижче цієї межі не опуститься жодна модель — навіть ідеальна")
розмір вибірки: 40 точок
шум: σ = 0.3, отже незвідна частина помилки σ² = 0.0900
нижче цієї межі не опуститься жодна модель — навіть ідеальна

1. Модель: поліном у базисі Лежандра¶

Складність моделі регулюємо степенем полінома. Але є технічна пастка, і про неї варто знати.

Якщо будувати матрицю ознак із сирих степенів $[1, x, x^2, \dots, x^{12}]$, її стовпці стають майже однаковими (усі виглядають як «щось, що росте»), матриця погано обумовлена, і на високих степенях результат псується чисельно, а не через перенавчання. Ми б міряли похибку арифметики замість дисперсії моделі.

Поліноми Лежандра дають той самий простір функцій, але їхні стовпці майже ортогональні. Простір той самий — отже, прогноз МНК має збігтися до останньої цифри. Це ми зараз і перевіримо.

In [2]:
def fit_polynomial(x_train, y_train, degree, x_new):
    """МНК-поліном заданого степеня в базисі Лежандра. Повертає прогноз у точках x_new."""
    design = legendre.legvander(x_train, degree)          # стовпці: P₀(x), P₁(x), …, P_d(x)
    coefficients, *_ = np.linalg.lstsq(design, y_train, rcond=None)
    return legendre.legvander(x_new, degree) @ coefficients


our_prediction = fit_polynomial(x_demo, y_demo, 7, x_test)
print(f"прогноз у перших трьох точках: {np.round(our_prediction[:3], 6)}")
прогноз у перших трьох точках: [1.522541 1.493871 1.462539]
In [3]:
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

# та сама задача, але через звичайні степені й sklearn
library_model = make_pipeline(PolynomialFeatures(7), LinearRegression())
library_model.fit(x_demo.reshape(-1, 1), y_demo)
library_prediction = library_model.predict(x_test.reshape(-1, 1))

print(f"найбільше розходження: {np.abs(our_prediction - library_prediction).max():.2e}")

assert np.allclose(our_prediction, library_prediction), "розрахунок розійшовся!"
print("\n✅ збігається — базис Лежандра дає той самий простір функцій,")
print("   але поводиться чисельно набагато краще на високих степенях")
найбільше розходження: 2.89e-15

✅ збігається — базис Лежандра дає той самий простір функцій,
   але поводиться чисельно набагато краще на високих степенях

2. Сорок паралельних світів¶

Ось головна ідея розділу 03 лекції: зміщення й дисперсія — властивості процедури навчання, а не однієї навченої моделі. Побачити їх на одному прогоні неможливо.

Тому влаштовуємо уявний експеримент по-справжньому: витягуємо багато різних навчальних вибірок з того самого джерела, на кожній навчаємо свою модель і дивимось на розкид прогнозів.

In [4]:
N_WORLDS = 300     # скільки паралельних світів проживаємо


def run_parallel_worlds(degree, size=SAMPLE_SIZE, worlds=N_WORLDS, seed=0):
    """Навчає `worlds` моделей на різних вибірках. Повертає матрицю прогнозів.

    Рядок — один світ, стовпець — одна тестова точка.
    """
    rng = np.random.default_rng(seed)
    predictions = np.zeros((worlds, len(x_test)))
    for world in range(worlds):
        x_train, y_train = draw_training_sample(rng, size)
        predictions[world] = fit_polynomial(x_train, y_train, degree, x_test)
    return predictions


predictions_degree_7 = run_parallel_worlds(degree=7)
print(f"матриця прогнозів: {predictions_degree_7.shape} (світів × тестових точок)")
print(f"у точці x = 0 прогнози гуляють від {predictions_degree_7[:, 85].min():.2f} "
      f"до {predictions_degree_7[:, 85].max():.2f}")
print(f"істина там: {f_test[85]:.2f}")
матриця прогнозів: (300, 171) (світів × тестових точок)
у точці x = 0 прогнози гуляють від -0.32 до 0.44
істина там: 0.00

Віяло¶

Кожна тонка лінія — модель з окремого світу. Товста рожева — середня модель $\bar{f}$, та сама, що стоїть у визначенні зміщення. Сірий пунктир — істина.

Дивись на дві речі окремо:

  • наскільки товста лінія розходиться з пунктиром — це зміщення;
  • наскільки широке віяло — це дисперсія.
In [5]:
fig, axes = plt.subplots(1, 3, figsize=(14, 4.2), sharey=True)

for ax, degree in zip(axes, [1, 5, 12]):
    predictions = run_parallel_worlds(degree)
    # малюємо лише 40 світів — інакше картинка перетвориться на суцільну пляму
    for world in range(40):
        ax.plot(x_test, predictions[world], color="teal", alpha=.16, lw=1)
    ax.plot(x_test, predictions.mean(axis=0), color="crimson", lw=3, label="середня модель")
    ax.plot(x_test, f_test, color="grey", ls="--", lw=2, label="істина")
    ax.set_title(f"степінь {degree}")
    ax.set_xlabel("x"); ax.grid(alpha=.25)

axes[0].set_ylabel("прогноз"); axes[0].set_ylim(-2.5, 2.5); axes[0].legend(loc="upper left")
plt.tight_layout(); plt.show()

print("Степінь 1: усі прямі лежать майже одна на одній (дисперсії немає),")
print("           але жодна навіть не намагається повторити форму істини.")
print("Степінь 12: середня модель лягла на істину чудово, зате окрема модель")
print("           може дати будь-що. Це і є компроміс, буквально очима.")
No description has been provided for this image
Степінь 1: усі прямі лежать майже одна на одній (дисперсії немає),
           але жодна навіть не намагається повторити форму істини.
Степінь 12: середня модель лягла на істину чудово, зате окрема модель
           може дати будь-що. Це і є компроміс, буквально очима.

3. Рахуємо три доданки числом¶

Формули з лекції, слово в слово:

$$\text{Зміщення}(x) = \mathbb{E}_D[\hat{f}_D(x)] - f(x), \qquad \text{Дисперсія}(x) = \mathbb{E}_D\big[(\hat{f}_D(x) - \mathbb{E}_D[\hat{f}_D(x)])^2\big]$$

Сподівання $\mathbb{E}_D$ береться по навчальних вибірках — тобто по рядках нашої матриці прогнозів. Потім усереднюємо по тестових точках, щоб отримати одне число на модель.

In [6]:
def decompose(predictions):
    """Розкладає помилку на три доданки. Повертає (зміщення², дисперсія, шум)."""
    average_model = predictions.mean(axis=0)                 # f̄(x): середнє по світах
    bias_squared = np.mean((average_model - f_test) ** 2)    # промах середньої моделі
    variance = np.mean(predictions.var(axis=0))              # розліт навколо середньої
    noise = NOISE_SIGMA ** 2                                 # підлога, однакова завжди
    return bias_squared, variance, noise


print(f"{'степінь':>8} {'зміщення²':>11} {'дисперсія':>11} {'шум':>8} {'разом':>9}")
for degree in [1, 5, 12]:
    bias_squared, variance, noise = decompose(run_parallel_worlds(degree))
    print(f"{degree:>8} {bias_squared:>11.4f} {variance:>11.4f} {noise:>8.4f} "
          f"{bias_squared + variance + noise:>9.4f}")
 степінь   зміщення²   дисперсія      шум     разом
       1      0.5298      0.0261   0.0900    0.6459
       5      0.0067      0.0131   0.0900    0.1098
      12      0.0012      0.5556   0.0900    0.6468

Перевірка 1: доданки справді сумуються¶

Лекція наполягає, що розклад — тотожність, а не наближення. Перевіримо це буквально: порахуємо середній квадрат відхилення прогнозів від істини напряму, без жодного розкладу, і порівняємо зі сумою зміщення² та дисперсії.

In [7]:
predictions = run_parallel_worlds(degree=7)
bias_squared, variance, noise = decompose(predictions)

# «в лоб»: середній по світах і по точках квадрат відхилення прогнозу від істини
straight_error = np.mean((predictions - f_test) ** 2)

print(f"зміщення² + дисперсія = {bias_squared + variance:.12f}")
print(f"пораховано напряму    = {straight_error:.12f}")

assert np.allclose(bias_squared + variance, straight_error), "розклад не сходиться!"
print("\n✅ тотожність підтверджена: доданків рівно два плюс шум, і вони не перекриваються")
зміщення² + дисперсія = 0.019430334275
пораховано напряму    = 0.019430334275

✅ тотожність підтверджена: доданків рівно два плюс шум, і вони не перекриваються

Перевірка 2: а тепер із справжнім шумом на тесті¶

Попередня перевірка була алгебраїчною: істину ми знали точно. Тепер зробимо чесніше — у кожному світі згенеруємо нові зашумлені тестові відповіді, як у реальному житті, і поміряємо помилку на них.

Тут уже працює статистика, тому точного збігу не буде — лише збіг у межах похибки Монте-Карло. Саме так і має бути.

In [8]:
noise_rng = np.random.default_rng(123)

# у кожному світі — свій свіжий шум на тестових точках
noisy_targets = f_test + noise_rng.normal(0, NOISE_SIGMA, predictions.shape)
measured_error = np.mean((noisy_targets - predictions) ** 2)

predicted_by_decomposition = bias_squared + variance + noise

print(f"поміряна помилка на зашумлених відповідях: {measured_error:.5f}")
print(f"передбачення розкладу (bias² + var + σ²):  {predicted_by_decomposition:.5f}")
print(f"розбіжність: {abs(measured_error - predicted_by_decomposition) / predicted_by_decomposition * 100:.2f}%")

assert np.allclose(measured_error, predicted_by_decomposition, rtol=0.05), \
    "розклад не описує реальну помилку!"
print("\n✅ розклад передбачає реальну помилку з точністю до похибки Монте-Карло")
поміряна помилка на зашумлених відповідях: 0.10888
передбачення розкладу (bias² + var + σ²):  0.10943
розбіжність: 0.50%

✅ розклад передбачає реальну помилку з точністю до похибки Монте-Карло

4. Розклад по складності — головна картинка теми¶

Ті самі числа, але для всіх степенів одразу. Сірий фундамент однаковий скрізь: це шум, підлога, нижче якої не опуститься ніхто.

In [9]:
degrees = np.arange(1, 13)
bias_by_degree = np.zeros(len(degrees))
variance_by_degree = np.zeros(len(degrees))

for i, degree in enumerate(degrees):
    bias_by_degree[i], variance_by_degree[i], _ = decompose(run_parallel_worlds(degree))

total_by_degree = bias_by_degree + variance_by_degree + NOISE_SIGMA ** 2
best_degree = degrees[total_by_degree.argmin()]

print(f"{'степінь':>8} {'зміщення²':>11} {'дисперсія':>11} {'разом':>9}")
for i, degree in enumerate(degrees):
    mark = "  ← мінімум" if degree == best_degree else ""
    print(f"{degree:>8} {bias_by_degree[i]:>11.4f} {variance_by_degree[i]:>11.4f} "
          f"{total_by_degree[i]:>9.4f}{mark}")
 степінь   зміщення²   дисперсія     разом
       1      0.5298      0.0261    0.6459
       2      0.5062      0.0421    0.6383
       3      0.1645      0.0267    0.2813
       4      0.1499      0.0449    0.2847
       5      0.0067      0.0131    0.1098
       6      0.0055      0.0177    0.1132
       7      0.0003      0.0192    0.1094  ← мінімум
       8      0.0003      0.0279    0.1182
       9      0.0002      0.0391    0.1292
      10      0.0003      0.1200    0.2104
      11      0.0005      0.2830    0.3735
      12      0.0012      0.5556    0.6468
In [10]:
fig, axes = plt.subplots(1, 2, figsize=(13, 4.4))

noise_row = np.full(len(degrees), NOISE_SIGMA ** 2)
axes[0].bar(degrees, noise_row, color="lightgrey", label="шум σ²")
axes[0].bar(degrees, bias_by_degree, bottom=noise_row, color="crimson", label="зміщення²")
axes[0].bar(degrees, variance_by_degree, bottom=noise_row + bias_by_degree,
            color="teal", label="дисперсія")
axes[0].set_yscale("log")
axes[0].set_title("З чого складається помилка (лог. шкала)")

axes[1].plot(degrees, bias_by_degree, "o-", color="crimson", lw=2, label="зміщення²")
axes[1].plot(degrees, variance_by_degree, "o-", color="teal", lw=2, label="дисперсія")
axes[1].plot(degrees, total_by_degree, "o-", color="black", lw=2.5, label="разом")
axes[1].axvline(best_degree, ls="--", color="grey", label=f"мінімум суми: {best_degree}")
axes[1].set_yscale("log")
axes[1].set_title("Зміщення падає, дисперсія росте")

for ax in axes:
    ax.set_xlabel("степінь полінома"); ax.legend(); ax.grid(alpha=.25)
plt.tight_layout(); plt.show()

print(f"Мінімум суми — на степені {best_degree}. Зверни увагу: це не там, де")
print("найменше зміщення, і не там, де найменша дисперсія. Оптимізувати треба суму.")
No description has been provided for this image
Мінімум суми — на степені 7. Зверни увагу: це не там, де
найменше зміщення, і не там, де найменша дисперсія. Оптимізувати треба суму.

5. Що бачить практик: криві train і test¶

Паралельних світів у житті немає — є одна вибірка. Що з неї видно?

Дві криві: помилка на навчальних даних і помилка на відкладених. Поводяться вони принципово по-різному, і саме розрив між ними ставить діагноз.

In [11]:
practice_rng = np.random.default_rng(2024)
x_train_one, y_train_one = draw_training_sample(practice_rng, size=SAMPLE_SIZE)

# відкладена вибірка з того самого джерела — те, що в житті називають тестом
x_holdout, y_holdout = draw_training_sample(practice_rng, size=400)

train_error = np.zeros(len(degrees))
holdout_error = np.zeros(len(degrees))

for i, degree in enumerate(degrees):
    prediction_on_train = fit_polynomial(x_train_one, y_train_one, degree, x_train_one)
    prediction_on_holdout = fit_polynomial(x_train_one, y_train_one, degree, x_holdout)
    train_error[i] = np.mean((y_train_one - prediction_on_train) ** 2)
    holdout_error[i] = np.mean((y_holdout - prediction_on_holdout) ** 2)

print(f"{'степінь':>8} {'train':>9} {'test':>9} {'розрив':>9}")
for i, degree in enumerate(degrees):
    print(f"{degree:>8} {train_error[i]:>9.4f} {holdout_error[i]:>9.4f} "
          f"{holdout_error[i] - train_error[i]:>9.4f}")
 степінь     train      test    розрив
       1    0.5289    0.5824    0.0535
       2    0.4833    0.5781    0.0948
       3    0.2317    0.2852    0.0535
       4    0.2310    0.2860    0.0550
       5    0.0799    0.1105    0.0306
       6    0.0794    0.1102    0.0308
       7    0.0669    0.0982    0.0313
       8    0.0646    0.0983    0.0337
       9    0.0426    0.1383    0.0958
      10    0.0422    0.1475    0.1054
      11    0.0378    0.1652    0.1274
      12    0.0376    0.1722    0.1346
In [12]:
fig, ax = plt.subplots(figsize=(9, 4.4))

ax.plot(degrees, train_error, "o-", color="crimson", lw=2.5, label="train")
ax.plot(degrees, holdout_error, "o-", color="teal", lw=2.5, label="test (відкладені)")
ax.axhline(NOISE_SIGMA ** 2, ls=":", color="grey", label="підлога σ² — нижче не буває")
ax.set_yscale("log")
ax.set_xlabel("степінь полінома"); ax.set_ylabel("MSE (лог. шкала)")
ax.set_title("Навчальна помилка падає завжди — тому вона нічого не каже")
ax.legend(); ax.grid(alpha=.25)
plt.tight_layout(); plt.show()

print("Діагностика за п'ять хвилин:")
print("  обидві криві високі й поруч → зміщення → ускладнювати")
print("  розрив великий              → дисперсія → спрощувати або добувати дані")
No description has been provided for this image
Діагностика за п'ять хвилин:
  обидві криві високі й поруч → зміщення → ускладнювати
  розрив великий              → дисперсія → спрощувати або добувати дані

6. Бутстреп: те саме, коли вибірка одна¶

Паралельні світи — уявний експеримент. Але дещо схоже можна зробити й насправді: бутстреп. Беремо одну наявну вибірку й багато разів витягуємо з неї $n$ обʼєктів з поверненням. Кожна така підвибірка трохи інша — от і «паралельні світи», зроблені з підручних матеріалів.

Бутстреп не знає істини, тому зміщення він оцінити не може. А от дисперсію — цілком: вона не потребує знання $f$.

In [13]:
def bootstrap_predictions(x_train, y_train, degree, resamples=N_WORLDS, seed=5):
    """Багато моделей, навчених на підвибірках з поверненням з однієї вибірки."""
    rng = np.random.default_rng(seed)
    size = len(x_train)
    predictions = np.zeros((resamples, len(x_test)))
    for i in range(resamples):
        # витягуємо номери рядків з поверненням — частина обʼєктів повториться,
        # частина не потрапить зовсім, і саме це створює різницю між моделями
        chosen = rng.integers(0, size, size)
        predictions[i] = fit_polynomial(x_train[chosen], y_train[chosen], degree, x_test)
    return predictions


print(f"{'степінь':>8} {'дисперсія (світи)':>19} {'дисперсія (бутстреп)':>22} {'у скільки разів':>17}")
for degree in [1, 5, 9]:
    _, ideal_variance, _ = decompose(run_parallel_worlds(degree))
    bootstrap_variance = np.mean(
        bootstrap_predictions(x_train_one, y_train_one, degree).var(axis=0))
    print(f"{degree:>8} {ideal_variance:>19.4f} {bootstrap_variance:>22.4f} "
          f"{bootstrap_variance / ideal_variance:>17.1f}")

print("\nНапрямок бутстреп ловить правильно: дисперсія росте зі складністю.")
print("А от абсолютні числа він завищує, і тим сильніше, чим складніша модель.")
print("Причина: підвибірка з поверненням містить у середньому лише 63% різних")
print("обʼєктів. Гнучкому поліному цього критично мало — він хапається за дублікати,")
print("і його розкид роздувається. Тому бутстреп — індикатор, а не вимірювальний прилад.")
 степінь   дисперсія (світи)   дисперсія (бутстреп)   у скільки разів
       1              0.0261                 0.0215               0.8
       5              0.0131                 0.0144               1.1
       9              0.0391                 0.2319               5.9

Напрямок бутстреп ловить правильно: дисперсія росте зі складністю.
А от абсолютні числа він завищує, і тим сильніше, чим складніша модель.
Причина: підвибірка з поверненням містить у середньому лише 63% різних
обʼєктів. Гнучкому поліному цього критично мало — він хапається за дублікати,
і його розкид роздувається. Тому бутстреп — індикатор, а не вимірювальний прилад.

7. Дані б'ють лише по дисперсії¶

Лекція: $\text{Дисперсія} \approx C/n$, а зміщення від $n$ не залежить взагалі. Перевіримо обидва твердження одразу, зафіксувавши складність.

In [14]:
sample_sizes = np.array([40, 80, 160, 320, 640])
bias_by_size = np.zeros(len(sample_sizes))
variance_by_size = np.zeros(len(sample_sizes))

for i, size in enumerate(sample_sizes):
    predictions = run_parallel_worlds(degree=7, size=size, worlds=200)
    bias_by_size[i], variance_by_size[i], _ = decompose(predictions)

print(f"{'n':>6} {'зміщення²':>11} {'дисперсія':>11} {'дисперсія × n':>15}")
for i, size in enumerate(sample_sizes):
    print(f"{size:>6} {bias_by_size[i]:>11.5f} {variance_by_size[i]:>11.5f} "
          f"{variance_by_size[i] * size:>15.3f}")

print("\nЗміщення² стоїть на місці — скільки даних не давай, поліном 7-го степеня")
print("лишиться поліномом 7-го степеня. Це і є сенс слова «систематична».")
print(f"\nА дисперсія × n майже не рухається: {variance_by_size[0] * sample_sizes[0]:.2f} → "
      f"{variance_by_size[-1] * sample_sizes[-1]:.2f}, тоді як сама n виросла в "
      f"{sample_sizes[-1] // sample_sizes[0]} разів.")
print("Тобто дисперсія падає приблизно як 1/n — трохи швидше, бо на малих вибірках")
print("до неї домішується чисельна нестійкість поліноміальної підгонки.")
     n   зміщення²   дисперсія   дисперсія × n
    40     0.00025     0.02039           0.816
    80     0.00014     0.00771           0.616
   160     0.00011     0.00358           0.573
   320     0.00012     0.00191           0.611
   640     0.00010     0.00086           0.549

Зміщення² стоїть на місці — скільки даних не давай, поліном 7-го степеня
лишиться поліномом 7-го степеня. Це і є сенс слова «систематична».

А дисперсія × n майже не рухається: 0.82 → 0.55, тоді як сама n виросла в 16 разів.
Тобто дисперсія падає приблизно як 1/n — трохи швидше, бо на малих вибірках
до неї домішується чисельна нестійкість поліноміальної підгонки.
In [15]:
fig, ax = plt.subplots(figsize=(8, 4.4))

ax.loglog(sample_sizes, variance_by_size, "o-", color="teal", lw=2.5, label="дисперсія")
ax.loglog(sample_sizes, bias_by_size + 1e-6, "o-", color="crimson", lw=2.5,
          label="зміщення² (+1e-6, щоб було видно на лог. шкалі)")
# еталонний нахил 1/n, прив'язаний до першої точки
reference = variance_by_size[0] * sample_sizes[0] / sample_sizes
ax.loglog(sample_sizes, reference, ls="--", color="grey", lw=2, label="еталон 1/n")

ax.set_xlabel("розмір навчальної вибірки n"); ax.set_ylabel("значення (лог. шкала)")
ax.set_title("Більше даних — менша дисперсія. Зміщення не рухається")
ax.legend(); ax.grid(alpha=.25, which="both")
plt.tight_layout(); plt.show()
No description has been provided for this image

8. Регуляризація: свідомий обмін¶

Найцікавіше. Візьмемо навмисно завелику модель — поліном 11-го степеня на 40 точках — і додамо штраф за величину коефіцієнтів:

$$L = \sum (y - \hat{y})^2 + \lambda \sum_{j \ge 1} \beta_j^2$$

Вільний член не штрафуємо: зсувати всю криву вгору-вниз — це не «зайва гнучкість», і карати за це немає за що.

Питання одне: чи виграш у дисперсії перевищить програш у зміщенні?

In [16]:
def run_worlds_with_ridge(lam, degree=11, worlds=N_WORLDS, seed=0):
    """Ті самі паралельні світи, але з L2-штрафом. Повертає матрицю прогнозів."""
    rng = np.random.default_rng(seed)
    penalty = np.eye(degree + 1)
    penalty[0, 0] = 0                      # вільний член не штрафуємо
    design_test = legendre.legvander(x_test, degree)

    predictions = np.zeros((worlds, len(x_test)))
    for world in range(worlds):
        x_train, y_train = draw_training_sample(rng)
        design = legendre.legvander(x_train, degree)
        # нормальне рівняння зі штрафом: (XᵀX + λI)β = Xᵀy
        coefficients = np.linalg.solve(design.T @ design + lam * penalty, design.T @ y_train)
        predictions[world] = design_test @ coefficients
    return predictions


lambdas = np.array([0.0, 1e-4, 1e-3, 1e-2, 1e-1, 1.0, 10.0])
bias_by_lambda = np.zeros(len(lambdas))
variance_by_lambda = np.zeros(len(lambdas))

for i, lam in enumerate(lambdas):
    bias_by_lambda[i], variance_by_lambda[i], _ = decompose(run_worlds_with_ridge(lam))

total_by_lambda = bias_by_lambda + variance_by_lambda + NOISE_SIGMA ** 2
best_lambda = lambdas[total_by_lambda.argmin()]

print(f"{'λ':>10} {'зміщення²':>11} {'дисперсія':>11} {'разом':>9}")
for i, lam in enumerate(lambdas):
    mark = "  ← найкраще" if lam == best_lambda else ""
    print(f"{lam:>10} {bias_by_lambda[i]:>11.4f} {variance_by_lambda[i]:>11.4f} "
          f"{total_by_lambda[i]:>9.4f}{mark}")

print(f"\nБез штрафу помилка {total_by_lambda[0]:.4f}, з найкращим λ — {total_by_lambda.min():.4f}.")
print(f"Виграш у {total_by_lambda[0] / total_by_lambda.min():.1f} раза — і це та сама модель,")
print("того самого 11-го степеня. Ми лише заборонили їй користуватися всією гнучкістю.")
         λ   зміщення²   дисперсія     разом
       0.0      0.0005      0.2830    0.3735
    0.0001      0.0001      0.0527    0.1428
     0.001      0.0001      0.0450    0.1351
      0.01      0.0001      0.0341    0.1241
       0.1      0.0005      0.0264    0.1169  ← найкраще
       1.0      0.0232      0.0188    0.1320
      10.0      0.2482      0.0152    0.3534

Без штрафу помилка 0.3735, з найкращим λ — 0.1169.
Виграш у 3.2 раза — і це та сама модель,
того самого 11-го степеня. Ми лише заборонили їй користуватися всією гнучкістю.
In [17]:
fig, ax = plt.subplots(figsize=(9, 4.4))

# нуль на логарифмічній осі не намалюєш, тому λ=0 показуємо окремою рискою
positive = lambdas > 0
ax.plot(lambdas[positive], bias_by_lambda[positive], "o-", color="crimson", lw=2.5,
        label="зміщення² — росте")
ax.plot(lambdas[positive], variance_by_lambda[positive], "o-", color="teal", lw=2.5,
        label="дисперсія — падає")
ax.plot(lambdas[positive], total_by_lambda[positive], "o-", color="black", lw=2.5,
        label="разом")
ax.axhline(total_by_lambda[0], ls=":", color="grey", label="разом при λ = 0")
ax.set_xscale("log"); ax.set_yscale("log")
ax.set_xlabel("λ (лог. шкала)"); ax.set_ylabel("внесок у помилку (лог. шкала)")
ax.set_title("λ обмінює дисперсію на зміщення — і сума має мінімум")
ax.legend(); ax.grid(alpha=.25, which="both")
plt.tight_layout(); plt.show()
No description has been provided for this image

9. Ті самі криві готовими інструментами¶

Усе, що ми рахували руками, у scikit-learn уже є двома функціями. Писати руками було потрібно, щоб зрозуміти, що саме вони рахують, — тепер можна користуватися.

  • validation_curve — помилка від складності (наша U-подібна крива);
  • learning_curve — помилка від розміру вибірки (та сама крива навчання, яку лекція радить будувати перед тим, як замовляти нові дані).

Обидві всередині роблять крос-валідацію, тому дають чесніші числа, ніж одне відкладене розбиття. Візьмемо трохи більшу вибірку — 150 точок, щоб кривій навчання було куди рости.

In [18]:
from sklearn.model_selection import validation_curve, learning_curve, KFold
from sklearn.preprocessing import StandardScaler

bigger_rng = np.random.default_rng(11)
x_big, y_big = draw_training_sample(bigger_rng, size=150)
x_big_column = x_big.reshape(-1, 1)          # sklearn чекає матрицю, а не вектор

# у пайплайні степінь стоїть окремим кроком — саме його ми й будемо крутити
polynomial_model = make_pipeline(PolynomialFeatures(), StandardScaler(), LinearRegression())
folds = KFold(n_splits=5, shuffle=True, random_state=0)

degree_range = np.arange(1, 15)
train_scores, test_scores = validation_curve(
    polynomial_model, x_big_column, y_big,
    param_name="polynomialfeatures__degree", param_range=degree_range,
    cv=folds, scoring="neg_mean_squared_error")

# sklearn повертає «чим більше, тим краще», тому міняємо знак назад на MSE
train_mse = -train_scores.mean(axis=1)
test_mse = -test_scores.mean(axis=1)
best_by_library = degree_range[test_mse.argmin()]

print(f"{'степінь':>8} {'train':>9} {'test':>9}")
for i, degree in enumerate(degree_range):
    mark = "  ← мінімум" if degree == best_by_library else ""
    print(f"{degree:>8} {train_mse[i]:>9.4f} {test_mse[i]:>9.4f}{mark}")

print(f"\nvalidation_curve на ОДНІЙ вибірці каже: степінь {best_by_library}.")
print(f"Наш чесний експеримент із {N_WORLDS} паралельними світами казав: степінь {best_degree}.")
print("Два зовсім різні шляхи привели до однієї відповіді — це і є найкраща перевірка.")
 степінь     train      test
       1    0.5721    0.5883
       2    0.5148    0.5400
       3    0.2591    0.2719
       4    0.2582    0.2763
       5    0.0813    0.0912
       6    0.0747    0.0834
       7    0.0654    0.0745  ← мінімум
       8    0.0650    0.0751
       9    0.0647    0.0752
      10    0.0645    0.0762
      11    0.0639    0.0775
      12    0.0614    0.0778
      13    0.0613    0.0799
      14    0.0606    0.0799

validation_curve на ОДНІЙ вибірці каже: степінь 7.
Наш чесний експеримент із 300 паралельними світами казав: степінь 7.
Два зовсім різні шляхи привели до однієї відповіді — це і є найкраща перевірка.
In [19]:
fig, axes = plt.subplots(1, 2, figsize=(13, 4.4))

axes[0].plot(degree_range, train_mse, "o-", color="crimson", lw=2.5, label="train")
axes[0].plot(degree_range, test_mse, "o-", color="teal", lw=2.5, label="test (CV)")
axes[0].axvline(best_by_library, ls="--", color="grey", label=f"мінімум: {best_by_library}")
axes[0].set_xlabel("степінь полінома"); axes[0].set_title("validation_curve: помилка від складності")

# крива навчання для замалої і для доречної моделі
for degree, colour in [(1, "crimson"), (5, "teal")]:
    sizes, _, holdout_scores = learning_curve(
        make_pipeline(PolynomialFeatures(degree), StandardScaler(), LinearRegression()),
        x_big_column, y_big, train_sizes=np.linspace(0.2, 1.0, 7),
        cv=folds, scoring="neg_mean_squared_error")
    axes[1].plot(sizes, -holdout_scores.mean(axis=1), "o-", color=colour, lw=2.5,
                 label=f"степінь {degree}")

axes[1].axhline(NOISE_SIGMA ** 2, ls=":", color="grey", label="підлога σ²")
axes[1].set_xlabel("розмір навчальної частини"); axes[1].set_title("learning_curve: помилка від кількості даних")

for ax in axes:
    ax.set_ylabel("MSE"); ax.legend(); ax.grid(alpha=.25)
plt.tight_layout(); plt.show()

print("Читаємо праву картинку так, як радить лекція:")
print("  степінь 1 — крива вийшла на плато високо над підлогою. Дані вже не")
print("              допоможуть, проблема в моделі. Замовляти розмітку — марно.")
print("  степінь 5 — крива ще спадає й тисне до підлоги σ². Ось тут нові дані")
print("              справді куплять тобі якість.")
No description has been provided for this image
Читаємо праву картинку так, як радить лекція:
  степінь 1 — крива вийшла на плато високо над підлогою. Дані вже не
              допоможуть, проблема в моделі. Замовляти розмітку — марно.
  степінь 5 — крива ще спадає й тисне до підлоги σ². Ось тут нові дані
              справді куплять тобі якість.

💻 Завдання¶

🟢 Рівень 1 — разом¶

  1. Зміни NOISE_SIGMA з 0.3 на 0.6 і перезапусти розділ 4. Який доданок змінився, а які лишились на місці? Чи зсунувся оптимальний степінь — і в який бік?
  2. Постав SAMPLE_SIZE = 15 і подивись на таблицю розкладу. На якому степені тепер починається катастрофа? Поясни зв'язок із кількістю параметрів моделі.

🟡 Рівень 2 — самостійно¶

  1. Заміни поліном на k-NN (sklearn.neighbors.KNeighborsRegressor) і побудуй ту саму таблицю зміщення²/дисперсії, але по $k$ від 1 до 25. Функція run_parallel_worlds майже не зміниться — треба лише інакше навчати модель.
  2. Побудуй криву навчання: помилка від розміру вибірки при фіксованій складності, для степенів 1 і 9 на одному графіку.

Зроблено, якщо: у k-NN зміщення й дисперсія поводяться дзеркально до полінома (мале $k$ = велика дисперсія), і ти написав(ла) одним реченням, чому саме так. А на кривій навчання видно, що для степеня 1 вона виходить на плато високо, а для степеня 9 продовжує падати.

🔴 Рівень 3 — виклик¶

  1. Візьми sklearn.ensemble.RandomForestRegressor і покажи числом те, про що говорить зірочка в таблиці розділу 10 лекції: при зростанні n_estimators від 1 до 200 дисперсія падає, а зміщення² майже не рухається.
  2. Порівняй це з GradientBoostingRegressor, у якого зростання n_estimators спершу знижує зміщення, а потім розганяє дисперсію.

Зроблено, якщо: є два графіки «зміщення² і дисперсія від кількості дерев» — для лісу й для бустингу — і по них видно, що ансамблі атакують різні доданки.