import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)

true_intercept = -4.3   # β₀ з лекції
true_slope = 0.36       # β₁ з лекції

n_students = 200
hours = rng.uniform(0, 40, n_students)

# істинна ймовірність скласти для кожного студента
probability_true = 1 / (1 + np.exp(-(true_intercept + true_slope * hours)))

# мітку кидаємо монеткою: так у даних зʼявляється чесна випадковість
passed = (rng.random(n_students) < probability_true).astype(int)

print(f"студентів: {n_students}, склали: {passed.sum()}, не склали: {n_students - passed.sum()}")
print(f"підготовка: від {hours.min():.1f} до {hours.max():.1f} годин")
print(f"істинна межа рішення: {-true_intercept / true_slope:.2f} години")

студентів: 200, склали: 139, не склали: 61
підготовка: від 0.3 до 39.7 годин
істинна межа рішення: 11.94 години

def sigmoid(z):
    """Стискає будь-яке число в проміжок (0, 1)."""
    # обрізаємо, щоб np.exp не переповнився на великих за модулем аргументах
    z = np.clip(z, -500, 500)
    return 1 / (1 + np.exp(-z))


z_grid = np.linspace(-8, 8, 401)

print(f"σ(0)          = {sigmoid(0):.6f}   (має бути 0.5)")
print(f"σ(-2)+σ(2)    = {sigmoid(-2) + sigmoid(2):.6f}   (має бути 1.0 — симетрія)")

# похідну рахуємо двома способами: за формулою і чисельно, різницею
derivative_formula = sigmoid(z_grid) * (1 - sigmoid(z_grid))
derivative_numeric = np.gradient(sigmoid(z_grid), z_grid)

print(f"максимум σ'   = {derivative_formula.max():.6f} при z = "
      f"{z_grid[derivative_formula.argmax()]:.2f}   (має бути 0.25 при z=0)")
print(f"формула vs чисельна похідна, найбільше розходження: "
      f"{np.abs(derivative_formula - derivative_numeric).max():.2e}")

σ(0)          = 0.500000   (має бути 0.5)
σ(-2)+σ(2)    = 1.000000   (має бути 1.0 — симетрія)
максимум σ'   = 0.250000 при z = 0.00   (має бути 0.25 при z=0)
формула vs чисельна похідна, найбільше розходження: 3.33e-05

fig, ax = plt.subplots(figsize=(9, 4))

ax.scatter(hours, passed, alpha=.45, s=40, color="teal", label="студенти (0 / 1)")
hours_line = np.linspace(0, 40, 300)
ax.plot(hours_line, sigmoid(true_intercept + true_slope * hours_line),
        color="crimson", lw=2.5, label="істинна ймовірність σ(z)")
ax.axvline(-true_intercept / true_slope, ls="--", color="grey", label="істинна межа")

ax.set_xlabel("годин підготовки"); ax.set_ylabel("склав (1) / не склав (0)")
ax.set_title("Дані та істинна залежність")
ax.legend(loc="center right"); ax.grid(alpha=.25)
plt.tight_layout(); plt.show()

# скільки коштує одна відповідь: даємо правильному класу різну впевненість
for probability_for_true_class in [0.99, 0.9, 0.5, 0.1, 0.01]:
    penalty = -np.log(probability_for_true_class)
    print(f"дав правильному класу {probability_for_true_class:>5} → штраф {penalty:6.3f}")

print("\nВисновок: впевнена помилка коштує необмежено дорого — саме тому")
print("модель, навчена на log-loss, воліє сказати «не знаю, 0.5», ніж збрехати.")

дав правильному класу  0.99 → штраф  0.010
дав правильному класу   0.9 → штраф  0.105
дав правильному класу   0.5 → штраф  0.693
дав правильному класу   0.1 → штраф  2.303
дав правильному класу  0.01 → штраф  4.605

Висновок: впевнена помилка коштує необмежено дорого — саме тому
модель, навчена на log-loss, воліє сказати «не знаю, 0.5», ніж збрехати.

def log_loss_by_hand(y_true, probabilities):
    """Бінарна крос-ентропія за означенням, усереднена по обʼєктах."""
    # обрізаємо: ln(0) = -inf зламав би всю подальшу арифметику
    p = np.clip(probabilities, 1e-15, 1 - 1e-15)
    return -np.mean(y_true * np.log(p) + (1 - y_true) * np.log(1 - p))


from sklearn.metrics import log_loss as sklearn_log_loss

# беремо істинні ймовірності як «прогноз» — аби було що міряти
our_value = log_loss_by_hand(passed, probability_true)
library_value = sklearn_log_loss(passed, probability_true)

print(f"наша       log-loss = {our_value:.10f}")
print(f"sklearn    log-loss = {library_value:.10f}")

assert np.allclose(our_value, library_value), "розрахунок розійшовся!"
print("\n✅ збігається")

наша       log-loss = 0.2464590070
sklearn    log-loss = 0.2464590070

✅ збігається

# ознаку стандартизуємо: інакше β₀ і β₁ живуть у різних масштабах
hours_scaled = (hours - hours.mean()) / hours.std()
design_matrix = np.c_[np.ones(n_students), hours_scaled]   # [1, x]


def log_loss_at(beta):
    """Log-loss як функція параметрів. Рахуємо через логіт — стійко й без обрізань."""
    logits = design_matrix @ beta
    # np.logaddexp(0, z) — це ln(1 + e^z), пораховане без переповнення
    return np.mean(np.logaddexp(0, logits) - passed * logits)


def mse_at(beta):
    """Та сама модель, але помилку міряємо квадратом відхилення від мітки."""
    return np.mean((sigmoid(design_matrix @ beta) - passed) ** 2)


print(f"log-loss у нулі: {log_loss_at(np.zeros(2)):.4f}   (ln 2 = {np.log(2):.4f})")
print(f"MSE      у нулі: {mse_at(np.zeros(2)):.4f}   (усім дали 0.5 → 0.25)")

log-loss у нулі: 0.6931   (ln 2 = 0.6931)
MSE      у нулі: 0.2500   (усім дали 0.5 → 0.25)

def how_much_above_chord(loss_function, beta_start, beta_end, steps=101):
    """Наскільки функція вилазить над хордою між двома наборами параметрів.

    Додатне число = порушення опуклості на цьому відрізку.
    """
    t_grid = np.linspace(0, 1, steps)
    values = np.zeros(steps)
    for i, t in enumerate(t_grid):
        values[i] = loss_function(beta_start + t * (beta_end - beta_start))
    chord = values[0] + t_grid * (values[-1] - values[0])
    return t_grid, values, chord


# проженемо багато випадкових напрямків і порахуємо, скільки з них порушують опуклість
chord_rng = np.random.default_rng(7)
n_directions = 300
violations_mse = 0
violations_log_loss = 0
worst_pair = None
worst_excess = 0.0

for _ in range(n_directions):
    start = chord_rng.uniform(-6, 6, 2)
    end = chord_rng.uniform(-6, 6, 2)

    _, values_mse, chord_mse = how_much_above_chord(mse_at, start, end)
    _, values_log, chord_log = how_much_above_chord(log_loss_at, start, end)

    excess_mse = np.max(values_mse - chord_mse)
    excess_log = np.max(values_log - chord_log)

    # 1e-9 — запас на похибку округлення, а не на математику
    if excess_mse > 1e-9:
        violations_mse += 1
    if excess_log > 1e-9:
        violations_log_loss += 1
    if excess_mse > worst_excess:
        worst_excess = excess_mse
        worst_pair = (start, end)

print(f"перевірено напрямків: {n_directions}")
print(f"MSE      вилізла над хордою у {violations_mse:>3} з {n_directions}")
print(f"log-loss вилізла над хордою у {violations_log_loss:>3} з {n_directions}")
print(f"\nнайбільше порушення в MSE: {worst_excess:.4f}")

перевірено напрямків: 300
MSE      вилізла над хордою у 206 з 300
log-loss вилізла над хордою у   0 з 300

найбільше порушення в MSE: 0.2009

start, end = worst_pair
t_grid, values_mse, chord_mse = how_much_above_chord(mse_at, start, end)
_, values_log, chord_log = how_much_above_chord(log_loss_at, start, end)

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

axes[0].plot(t_grid, values_mse, color="crimson", lw=2.5, label="MSE вздовж відрізка")
axes[0].plot(t_grid, chord_mse, ls="--", color="grey", lw=2, label="хорда")
axes[0].fill_between(t_grid, chord_mse, values_mse,
                     where=values_mse > chord_mse, color="crimson", alpha=.18)
axes[0].set_title("MSE: крива вилазить над хордою → невипукла")

axes[1].plot(t_grid, values_log, color="teal", lw=2.5, label="log-loss вздовж відрізка")
axes[1].plot(t_grid, chord_log, ls="--", color="grey", lw=2, label="хорда")
axes[1].set_title("log-loss: ніде не вище за хорду → випукла")

for ax in axes:
    ax.set_xlabel("положення на відрізку між двома наборами параметрів")
    ax.legend(); ax.grid(alpha=.25)
plt.tight_layout(); plt.show()

print("Практичний наслідок: у MSE спуск може застрягти в плато або локальній ямі,")
print("у log-loss мінімум рівно один — куди б ти не стартував, прийдеш туди ж.")

Практичний наслідок: у MSE спуск може застрягти в плато або локальній ямі,
у log-loss мінімум рівно один — куди б ти не стартував, прийдеш туди ж.

def log_loss_gradient(beta):
    """Градієнт log-loss за формулою з лекції."""
    predicted = sigmoid(design_matrix @ beta)
    return design_matrix.T @ (predicted - passed) / n_students


def numeric_gradient(loss_function, beta, step=1e-6):
    """Похідна «в лоб»: зсуваємо кожен параметр окремо й дивимось на приріст втрати."""
    gradient = np.zeros_like(beta)
    for i in range(len(beta)):
        shifted_up = beta.copy()
        shifted_down = beta.copy()
        shifted_up[i] += step
        shifted_down[i] -= step
        gradient[i] = (loss_function(shifted_up) - loss_function(shifted_down)) / (2 * step)
    return gradient


test_point = np.array([0.7, -1.3])
by_formula = log_loss_gradient(test_point)
by_numbers = numeric_gradient(log_loss_at, test_point)

print(f"градієнт за формулою:  {by_formula}")
print(f"градієнт чисельно:     {by_numbers}")

assert np.allclose(by_formula, by_numbers, atol=1e-7), "формула градієнта неправильна!"
print("\n✅ формула Xᵀ(σ(Xβ) − y) підтверджена чисельно")

градієнт за формулою:  [-0.0710591  -0.58271979]
градієнт чисельно:     [-0.0710591  -0.58271979]

✅ формула Xᵀ(σ(Xβ) − y) підтверджена чисельно

def fit_logistic_by_gradient_descent(eta=0.5, steps=20000):
    """Навчання логістичної регресії з нуля. Повертає коефіцієнти та історію втрат."""
    beta = np.zeros(2)
    loss_history = np.zeros(steps)
    for step in range(steps):
        loss_history[step] = log_loss_at(beta)
        beta = beta - eta * log_loss_gradient(beta)
    return beta, loss_history


our_beta, loss_history = fit_logistic_by_gradient_descent()

print(f"β₀ = {our_beta[0]:.6f}")
print(f"β₁ = {our_beta[1]:.6f}   (у стандартизованих одиницях)")
print(f"log-loss на старті: {loss_history[0]:.6f}")
print(f"log-loss у кінці:   {loss_history[-1]:.6f}")

β₀ = 2.542088
β₁ = 3.697196   (у стандартизованих одиницях)
log-loss на старті: 0.693147
log-loss у кінці:   0.245480

fig, ax = plt.subplots(figsize=(9, 4))

for eta in [0.01, 0.1, 0.5, 2.0]:
    _, history = fit_logistic_by_gradient_descent(eta=eta, steps=300)
    ax.plot(history, lw=2, label=f"η = {eta}")

ax.axhline(loss_history[-1], color="grey", ls="--", lw=1.5, label="досягнутий мінімум")
ax.set_xlabel("ітерація"); ax.set_ylabel("log-loss")
ax.set_title("Швидкість збіжності залежно від кроку")
ax.legend(); ax.grid(alpha=.25)
plt.tight_layout(); plt.show()

from sklearn.linear_model import LogisticRegression

# penalty=None — вимикаємо регуляризацію, інакше коефіцієнти навмисно не збіжаться
library_model = LogisticRegression(penalty=None, solver="lbfgs", max_iter=10000, tol=1e-10)
library_model.fit(hours_scaled.reshape(-1, 1), passed)

library_beta = np.array([library_model.intercept_[0], library_model.coef_[0, 0]])

print(f"наш спуск : β₀={our_beta[0]:.8f}  β₁={our_beta[1]:.8f}")
print(f"sklearn   : β₀={library_beta[0]:.8f}  β₁={library_beta[1]:.8f}")
print(f"різниця   : {np.abs(our_beta - library_beta)}")

assert np.allclose(our_beta, library_beta, atol=1e-6), "коефіцієнти розійшлися!"
print("\n✅ збігається — всередині бібліотеки той самий градієнтний спуск,")
print("   тільки розумніший (lbfgs збігається за десятки ітерацій замість тисяч)")

наш спуск : β₀=2.54208823  β₁=3.69719577
sklearn   : β₀=2.54208823  β₁=3.69719577
різниця   : [1.25742439e-09 5.82824011e-10]

✅ збігається — всередині бібліотеки той самий градієнтний спуск,
   тільки розумніший (lbfgs збігається за десятки ітерацій замість тисяч)

mean_hours = hours.mean()
std_hours = hours.std()

slope_in_hours = our_beta[1] / std_hours
intercept_in_hours = our_beta[0] - our_beta[1] * mean_hours / std_hours

print(f"відновили: β₀ = {intercept_in_hours:7.3f}   β₁ = {slope_in_hours:.3f}")
print(f"істина   : β₀ = {true_intercept:7.3f}   β₁ = {true_slope:.3f}")
print(f"\nмежа рішення: наша {-intercept_in_hours / slope_in_hours:.2f} год, "
      f"істинна {-true_intercept / true_slope:.2f} год")
print("\nРозбіжність — не помилка коду, а ціна скінченної вибірки:")
print("200 студентів дають лише наближення до істини.")

відновили: β₀ =  -3.942   β₁ = 0.328
істина   : β₀ =  -4.300   β₁ = 0.360

межа рішення: наша 12.01 год, істинна 11.94 год

Розбіжність — не помилка коду, а ціна скінченної вибірки:
200 студентів дають лише наближення до істини.

predicted_probability = sigmoid(design_matrix @ our_beta)


def confusion_metrics(y_true, probabilities, threshold):
    """Precision, recall і accuracy при заданому порозі."""
    predicted_label = (probabilities >= threshold).astype(int)
    true_positive = np.sum((predicted_label == 1) & (y_true == 1))
    false_positive = np.sum((predicted_label == 1) & (y_true == 0))
    false_negative = np.sum((predicted_label == 0) & (y_true == 1))
    # якщо модель нікого не позначила позитивним, precision не визначена — беремо 0
    if true_positive + false_positive == 0:
        precision = 0.0
    else:
        precision = true_positive / (true_positive + false_positive)

    recall = true_positive / (true_positive + false_negative)
    accuracy = np.mean(predicted_label == y_true)
    return precision, recall, accuracy


print(f"{'поріг':>7} {'precision':>10} {'recall':>8} {'accuracy':>9}")
for threshold in [0.1, 0.3, 0.5, 0.7, 0.9]:
    precision, recall, accuracy = confusion_metrics(passed, predicted_probability, threshold)
    print(f"{threshold:>7.2f} {precision:>10.3f} {recall:>8.3f} {accuracy:>9.3f}")

  поріг  precision   recall  accuracy
   0.10      0.793    0.993     0.815
   0.30      0.881    0.957     0.880
   0.50      0.921    0.928     0.895
   0.70      0.961    0.878     0.890
   0.90      0.981    0.734     0.805

from sklearn.metrics import precision_score, recall_score

our_precision, our_recall, _ = confusion_metrics(passed, predicted_probability, 0.5)
label_at_half = (predicted_probability >= 0.5).astype(int)

library_precision = precision_score(passed, label_at_half)
library_recall = recall_score(passed, label_at_half)

print(f"наші    : precision={our_precision:.6f}  recall={our_recall:.6f}")
print(f"sklearn : precision={library_precision:.6f}  recall={library_recall:.6f}")

assert np.allclose([our_precision, our_recall], [library_precision, library_recall]), \
    "метрики розійшлися!"
print("\n✅ збігається")

наші    : precision=0.921429  recall=0.928058
sklearn : precision=0.921429  recall=0.928058

✅ збігається

thresholds = np.linspace(0.01, 0.99, 199)
precision_curve = np.zeros_like(thresholds)
recall_curve = np.zeros_like(thresholds)

for i, threshold in enumerate(thresholds):
    precision_curve[i], recall_curve[i], _ = confusion_metrics(
        passed, predicted_probability, threshold)

required_recall = 0.90
suitable = thresholds[recall_curve >= required_recall]
chosen_threshold = suitable.max()
chosen_precision, chosen_recall, _ = confusion_metrics(
    passed, predicted_probability, chosen_threshold)

print(f"вимога: recall ≥ {required_recall}")
print(f"обраний поріг: {chosen_threshold:.3f}")
print(f"при ньому recall = {chosen_recall:.3f}, precision = {chosen_precision:.3f}")

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(thresholds, precision_curve, color="crimson", lw=2.5, label="precision")
ax.plot(thresholds, recall_curve, color="teal", lw=2.5, label="recall")
ax.axhline(required_recall, ls=":", color="grey", label=f"вимога recall ≥ {required_recall}")
ax.axvline(chosen_threshold, ls="--", color="black", lw=1.5, label="обраний поріг")
ax.set_xlabel("поріг"); ax.set_ylabel("значення метрики")
ax.set_title("Поріг торгує precision на recall")
ax.legend(); ax.grid(alpha=.25)
plt.tight_layout(); plt.show()

вимога: recall ≥ 0.9
обраний поріг: 0.540
при ньому recall = 0.906, precision = 0.926

odds_ratio = np.exp(slope_in_hours)
print(f"β₁ = {slope_in_hours:.4f} → e^β₁ = {odds_ratio:.4f}")
print(f"Кожна додаткова година множить шанси скласти на {odds_ratio:.3f}, "
      f"тобто підвищує їх на {(odds_ratio - 1) * 100:.1f}%.\n")

checkpoints = np.array([0, 5, 10, 15, 20, 25, 30])
probability_at = sigmoid(intercept_in_hours + slope_in_hours * checkpoints)
odds_at = probability_at / (1 - probability_at)

print(f"{'годин':>6} {'p':>8} {'шанси':>10} {'×шанси':>9} {'Δp':>8}")
for i in range(len(checkpoints)):
    if i == 0:
        print(f"{checkpoints[i]:>6} {probability_at[i]:>8.4f} {odds_at[i]:>10.4f} "
              f"{'—':>9} {'—':>8}")
    else:
        print(f"{checkpoints[i]:>6} {probability_at[i]:>8.4f} {odds_at[i]:>10.4f} "
              f"{odds_at[i] / odds_at[i - 1]:>9.4f} "
              f"{probability_at[i] - probability_at[i - 1]:>8.4f}")

β₁ = 0.3282 → e^β₁ = 1.3885
Кожна додаткова година множить шанси скласти на 1.389, тобто підвищує їх на 38.9%.

 годин        p      шанси    ×шанси       Δp
     0   0.0190     0.0194         —        —
     5   0.0911     0.1002    5.1610   0.0720
    10   0.3408     0.5170    5.1610   0.2498
    15   0.7274     2.6683    5.1610   0.3866
    20   0.9323    13.7711    5.1610   0.2049
    25   0.9861    71.0726    5.1610   0.0538
    30   0.9973   366.8046    5.1610   0.0112

# множник на шанси має бути сталим, а приріст імовірності — ні
odds_multipliers = odds_at[1:] / odds_at[:-1]
probability_steps = np.diff(probability_at)

expected_multiplier = np.exp(slope_in_hours * 5)   # крок таблиці — 5 годин

print(f"множники на шанси: {np.round(odds_multipliers, 6)}")
print(f"очікуваний e^(5β₁) = {expected_multiplier:.6f}")

assert np.allclose(odds_multipliers, expected_multiplier), "множник шансів не сталий!"
print("\n✅ множник на шанси сталий на кожному кроці")

print(f"\nа приріст імовірності гуляє: від {probability_steps.min():.4f} "
      f"до {probability_steps.max():.4f} — різниця у "
      f"{probability_steps.max() / probability_steps.min():.0f} разів")
print("Тому коефіцієнт логістичної регресії читають лише через шанси.")

множники на шанси: [5.160983 5.160983 5.160983 5.160983 5.160983 5.160983]
очікуваний e^(5β₁) = 5.160983

✅ множник на шанси сталий на кожному кроці

а приріст імовірності гуляє: від 0.0112 до 0.3866 — різниця у 35 разів
Тому коефіцієнт логістичної регресії читають лише через шанси.

Практика 09 · Логістична регресія¶

0. Дані¶

1. Сигмоїда руками¶

Як виглядають дані й істинна крива¶

2. Log-loss руками¶

3. Чому не MSE: тест хорди¶

4. Градієнтний спуск з нуля¶

Крок η вирішує все¶

Звірка зі scikit-learn¶

Повертаємось у години¶

5. Поріг — важіль поза моделлю¶

Як обирають поріг у житті¶

6. Коефіцієнт як відношення шансів¶

💻 Завдання¶

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

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

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