Практика 13 · Зміщення проти дисперсії¶
📖 Теорія: відкрий
lecture.htmlу цій же теці. 📝 Домашнє завдання:homework.md. 🧠 Тест:quiz.html.
Лекція каже: помилка розкладається на три доданки — шум, зміщення² і дисперсію, і цей розклад точний. Тут ми його порахуємо руками й переконаємось, що доданки справді сумуються в загальну помилку.
Що зробимо:
- Створимо «сорок паралельних світів» — багато навчальних вибірок з одного джерела
- Побудуємо віяло прогнозів і побачимо дисперсію очима
- Порахуємо bias², variance і шум числом та перевіримо, що вони складаються в очікувану помилку
- Побудуємо криві train/test від складності — те, що бачить практик з однією вибіркою
- Зробимо те саме через бутстреп, коли вибірка одна-єдина
- Покажемо, що дані б'ють лише по дисперсії, і то як $1/n$
- Обміняємо дисперсію на зміщення регуляризацією — і виграємо в сумі
- Побудуємо те саме готовими інструментами
learning_curveіvalidation_curve
0. Істина, шум і вибірка¶
Влаштовуємо світ, у якому ми знаємо істину. У житті так не буває — саме тому bias і variance неможливо поміряти на реальних даних. Тут можна, бо істину задали ми самі.
$$y = f(x) + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0,\ \sigma^2)$$
Функція $f$ навмисно не поліном: жоден скінченний степінь не відтворить її точно, тому зміщення буде справжнім, а не штучним.
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}]$, її стовпці стають майже однаковими (усі виглядають як «щось, що росте»), матриця погано обумовлена, і на високих степенях результат псується чисельно, а не через перенавчання. Ми б міряли похибку арифметики замість дисперсії моделі.
Поліноми Лежандра дають той самий простір функцій, але їхні стовпці майже ортогональні. Простір той самий — отже, прогноз МНК має збігтися до останньої цифри. Це ми зараз і перевіримо.
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]
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 лекції: зміщення й дисперсія — властивості процедури навчання, а не однієї навченої моделі. Побачити їх на одному прогоні неможливо.
Тому влаштовуємо уявний експеримент по-справжньому: витягуємо багато різних навчальних вибірок з того самого джерела, на кожній навчаємо свою модель і дивимось на розкид прогнозів.
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}$, та сама, що стоїть у визначенні зміщення. Сірий пунктир — істина.
Дивись на дві речі окремо:
- наскільки товста лінія розходиться з пунктиром — це зміщення;
- наскільки широке віяло — це дисперсія.
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(" може дати будь-що. Це і є компроміс, буквально очима.")
Степінь 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$ береться по навчальних вибірках — тобто по рядках нашої матриці прогнозів. Потім усереднюємо по тестових точках, щоб отримати одне число на модель.
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: доданки справді сумуються¶
Лекція наполягає, що розклад — тотожність, а не наближення. Перевіримо це буквально: порахуємо середній квадрат відхилення прогнозів від істини напряму, без жодного розкладу, і порівняємо зі сумою зміщення² та дисперсії.
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: а тепер із справжнім шумом на тесті¶
Попередня перевірка була алгебраїчною: істину ми знали точно. Тепер зробимо чесніше — у кожному світі згенеруємо нові зашумлені тестові відповіді, як у реальному житті, і поміряємо помилку на них.
Тут уже працює статистика, тому точного збігу не буде — лише збіг у межах похибки Монте-Карло. Саме так і має бути.
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. Розклад по складності — головна картинка теми¶
Ті самі числа, але для всіх степенів одразу. Сірий фундамент однаковий скрізь: це шум, підлога, нижче якої не опуститься ніхто.
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
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("найменше зміщення, і не там, де найменша дисперсія. Оптимізувати треба суму.")
Мінімум суми — на степені 7. Зверни увагу: це не там, де найменше зміщення, і не там, де найменша дисперсія. Оптимізувати треба суму.
5. Що бачить практик: криві train і test¶
Паралельних світів у житті немає — є одна вибірка. Що з неї видно?
Дві криві: помилка на навчальних даних і помилка на відкладених. Поводяться вони принципово по-різному, і саме розрив між ними ставить діагноз.
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
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(" розрив великий → дисперсія → спрощувати або добувати дані")
Діагностика за п'ять хвилин: обидві криві високі й поруч → зміщення → ускладнювати розрив великий → дисперсія → спрощувати або добувати дані
6. Бутстреп: те саме, коли вибірка одна¶
Паралельні світи — уявний експеримент. Але дещо схоже можна зробити й насправді: бутстреп. Беремо одну наявну вибірку й багато разів витягуємо з неї $n$ обʼєктів з поверненням. Кожна така підвибірка трохи інша — от і «паралельні світи», зроблені з підручних матеріалів.
Бутстреп не знає істини, тому зміщення він оцінити не може. А от дисперсію — цілком: вона не потребує знання $f$.
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$ не залежить взагалі. Перевіримо обидва твердження одразу, зафіксувавши складність.
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 — трохи швидше, бо на малих вибірках
до неї домішується чисельна нестійкість поліноміальної підгонки.
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()
8. Регуляризація: свідомий обмін¶
Найцікавіше. Візьмемо навмисно завелику модель — поліном 11-го степеня на 40 точках — і додамо штраф за величину коефіцієнтів:
$$L = \sum (y - \hat{y})^2 + \lambda \sum_{j \ge 1} \beta_j^2$$
Вільний член не штрафуємо: зсувати всю криву вгору-вниз — це не «зайва гнучкість», і карати за це немає за що.
Питання одне: чи виграш у дисперсії перевищить програш у зміщенні?
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-го степеня. Ми лише заборонили їй користуватися всією гнучкістю.
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()
9. Ті самі криві готовими інструментами¶
Усе, що ми рахували руками, у scikit-learn уже є двома функціями. Писати руками
було потрібно, щоб зрозуміти, що саме вони рахують, — тепер можна користуватися.
validation_curve— помилка від складності (наша U-подібна крива);learning_curve— помилка від розміру вибірки (та сама крива навчання, яку лекція радить будувати перед тим, як замовляти нові дані).
Обидві всередині роблять крос-валідацію, тому дають чесніші числа, ніж одне відкладене розбиття. Візьмемо трохи більшу вибірку — 150 точок, щоб кривій навчання було куди рости.
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.
Два зовсім різні шляхи привели до однієї відповіді — це і є найкраща перевірка.
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(" справді куплять тобі якість.")
Читаємо праву картинку так, як радить лекція:
степінь 1 — крива вийшла на плато високо над підлогою. Дані вже не
допоможуть, проблема в моделі. Замовляти розмітку — марно.
степінь 5 — крива ще спадає й тисне до підлоги σ². Ось тут нові дані
справді куплять тобі якість.
💻 Завдання¶
🟢 Рівень 1 — разом¶
- Зміни
NOISE_SIGMAз0.3на0.6і перезапусти розділ 4. Який доданок змінився, а які лишились на місці? Чи зсунувся оптимальний степінь — і в який бік? - Постав
SAMPLE_SIZE = 15і подивись на таблицю розкладу. На якому степені тепер починається катастрофа? Поясни зв'язок із кількістю параметрів моделі.
🟡 Рівень 2 — самостійно¶
- Заміни поліном на k-NN (
sklearn.neighbors.KNeighborsRegressor) і побудуй ту саму таблицю зміщення²/дисперсії, але по $k$ від 1 до 25. Функціяrun_parallel_worldsмайже не зміниться — треба лише інакше навчати модель. - Побудуй криву навчання: помилка від розміру вибірки при фіксованій складності, для степенів 1 і 9 на одному графіку.
Зроблено, якщо: у k-NN зміщення й дисперсія поводяться дзеркально до полінома (мале $k$ = велика дисперсія), і ти написав(ла) одним реченням, чому саме так. А на кривій навчання видно, що для степеня 1 вона виходить на плато високо, а для степеня 9 продовжує падати.
🔴 Рівень 3 — виклик¶
- Візьми
sklearn.ensemble.RandomForestRegressorі покажи числом те, про що говорить зірочка в таблиці розділу 10 лекції: при зростанніn_estimatorsвід 1 до 200 дисперсія падає, а зміщення² майже не рухається. - Порівняй це з
GradientBoostingRegressor, у якого зростанняn_estimatorsспершу знижує зміщення, а потім розганяє дисперсію.
Зроблено, якщо: є два графіки «зміщення² і дисперсія від кількості дерев» — для лісу й для бустингу — і по них видно, що ансамблі атакують різні доданки.