Семестр 3 · Модуль 3 · Неделя 8

SciPy: численные расчёты и честная визуализация

quad · odeint / solve_ivp · minimize · curve_fit · interp1d · scipy.stats · принципы честной визуализации

2 ч теории. Данные — reaction.csv (25 экспериментальных точек). После лекции — live-coding и задачи семинара.

Прежде чем начать: как бы вы «проверили закон» по точкам с шумом — например, зависимость скорости реакции от температуры?

Подбор кривой — как физики проверяют законы

По экспериментальным точкам с шумом восстановить параметры модели — метод Методик вычислений и Исследования операций.
Так Аррениус проверял зависимость скорости реакции от температуры (экспонента); так подбирают затухание маятника, рост популяции.
Интегралы и дифференциальные уравнения — база следующих курсов численных методов.

reaction.csv: 25 точек (x 0..10 с, y 5.5..0.7 отн. ед.) и экспонента A·exp(−λt) + c — то, что восстановим на слайде 13.

Что за «модель» вы бы выбрали для данных, которые сначала быстро падают, потом выходят на плато? (экспонента + константа)

Retrieval по неделе 7: savefig и подписи

Не открывая конспект: в каком порядке относительно plt.show() вызывается fig.savefig, иначе файл выйдет пустым? (до)
Перечислите минимум подписей, без которых график «не считается сданным». (title, оси с единицами, легенда при >1 серии)
Ось Y bar-диаграммы начинается с… (нуля)

Это не повторение ради повторения — сегодня мы продолжим тему честных графиков и доделаем её полностью.

Где должен стоять savefig относительно show, чтобы файл не вышел пустым?

К концу лекции вы сможете…

Вычислять интегралы через quad и понимать оценку ошибки;
Решать ОДУ распада и роста через odeint;
Подбирать параметры модели к данным через curve_fit (p0, pcov, RMSE, остатки);
Отличать интерполяцию от экстраполяции и объяснять, почему экстраполяция опасна;
Применять scipy.stats и честно трактовать p-value;
Узнавать и «чинить» манипулирующие графики.

Где в этих целях «ядро» блока 2 вашего проекта модуля? (curve_fit)

quad возвращает (значение, ошибка): ∫₀¹ x² dx = 1/3

from scipy.integrate import quad
quad(lambda x: x**2, 0, 1)
(0.33333333333333337, 3.7e-15)

Правило: quad(f, a, b) возвращает кортеж (значение, оценка ошибки). ∫₀¹ x² dx = 1/3 — сверили с аналитикой: численный результат совпал с точностью до 10⁻¹⁵. Всегда показывайте оба числа: без оценки ошибки значение неполно.

Как проверить, что quad посчитал правильно, если аналитического ответа нет? (известные интегралы для самопроверки; другой метод; оценка ошибки из кортежа)

Площадь под нормальной кривой ≈ 1 — свойство плотности

from scipy.stats import norm
quad(norm.pdf, -np.inf, np.inf)
(0.9999999999999998, 1.0e-08)

Правило: площадь под функцией плотности вероятности равна 1 — это её определяющее свойство; norm.pdf — плотность стандартного нормального распределения. Не путать: norm.pdf (плотность) ≠ norm.cdf (функция распределения).

Вы не знаете аналитику — но знаете свойство: результат обязан быть 1. Что говорит то, что quad вернул 1 с ошибкой 1e-8?

Площадь под кривой затухания = 20.84 «отн. ед. × с»

def model(t, A, lam, c):
    return A * np.exp(-lam * t) + c

area, _ = quad(lambda t: model(t, 5.0, 0.3, 0.5), 0, 10)
print(area)   # 20.837
20.837

Смысл: площадь под кривой затухания от 0 до 10 с — 20.84 «отн. ед. × с» — например, суммарное количество реагирующего вещества за первые 10 секунд.

Здесь параметры заданы (A=5, λ=0.3, c=0.5) «вручную». Как мы получим их из данных дальше в лекции?

odeint: экспоненциальный распад y′ = −k·y, сверка с аналитикой

from scipy.integrate import odeint

def decay(y, t, k):
    return -k * y          # y' = -k·y

t = np.linspace(0, 10, 50)
sol = odeint(decay, 5.0, t, args=(0.3,))
sol[[0, 25, -1], 0]
array([5.   , 1.082, 0.249])

Правило: функция правой части принимает (y, t, …) и возвращает производную; odeint(f, y0, t, args=…) интегрирует от начального условия y0. Сверка: y(10) = 5·exp(−0.3·10) ≈ 0.249 — совпало. Порядок аргументов f(y, t), а не f(t, y)!

Что означает выборка sol[[0, 25, -1], 0]? (анонс predict)

odeint: логистический рост → плато K

def logistic(y, t, r, K):
    return r * y * (1 - y / K)

t2 = np.linspace(0, 20, 60)
sol2 = odeint(logistic, 0.1, t2, args=(0.5, 10.0))
sol2[[0, 30, -1], 0]
array([0.1 , 6.2 , 9.955])

Правило: логистический рост y′ = r·y·(1−y/K): при малых y растёт почти экспоненциально, при приближении к вместимости K тормозится и выходит на плато K=10.

Модель, не данные: распад y(0)=5, k=0.3; логистический рост r=0.5, K=10 — плато K=10.

Что, если K = 10, а y0 = 20? Куда пойдёт кривая?

SIR-модель — зачем решать ОДУ на самом деле

SIR-модель эпидемии: три уравнения для Susceptible / Infected / Recovered.
S′ = −β·S·I, I′ = β·S·I − γ·I, R′ = γ·I — те же приёмы: odeint с системой из 3 уравнений.
Такие модели считают в эпидемиологии, экологии, экономике.

Сегодня мы решали одно уравнение — система решается ровно так же: функция правой части возвращает вектор производных.

Сколько начальных условий нужно системе из трёх уравнений? (три)

Predict: sol[[0, 25, -1], 0]

t = np.linspace(0, 10, 50)
sol = odeint(decay, 5.0, t, args=(0.3,))
print(sol.shape)
print(sol[[0, 25, -1], 0])

Подсказка: odeint поддерживает системы — подумайте, какой размерности будет результат.

Не запуская: какая форма у sol? Что выведет вторая строка?

minimize на функции Розенброка: минимум (1, 1)

from scipy.optimize import minimize, rosen
minimize(rosen, x0=[0, 0], method='BFGS')
fun: 2.84e-11
x:   [1.0  1.0]

Правило: minimize ищет минимум многомерной функции; функция Розенброка имеет минимум в точке (1, 1) — BFGS нашёл его с точностью до 10⁻¹¹. Методы: Nelder-Mead (без производных), BFGS (градиентные). curve_fit внутри делает ровно это: минимизирует сумму квадратов остатков.

Что произойдёт, если стартовать minimize из далёкой точки? (может сойтись в другой локальный минимум — для Розенброка риск велик)

curve_fit к reaction.csv: A=4.94±0.10, λ=0.30±0.02, c=0.53±0.09

from scipy.optimize import curve_fit

def model(t, A, lam, c):
    return A * np.exp(-lam * t) + c

popt, pcov = curve_fit(model, x, y, p0=[5.0, 0.3, 0.5])
perr = np.sqrt(np.diag(pcov))
resid = y - model(x, *popt)
rmse = np.sqrt(np.mean(resid**2))
A = 4.94 ± 0.10
λ = 0.30 ± 0.02
c = 0.53 ± 0.09
RMSE = 0.121

Правило: curve_fit — метод наименьших квадратов; возвращает (popt, pcov), где pcov — ковариационная матрица; perr = sqrt(diag(pcov)) — ошибки параметров. Истинные (5, 0.3, 0.5) внутри ±2σ. RMSE=0.121 мал относительно размаха данных (5.5..0.7).

Почему параметры всегда показывают с ошибками (±)? Что изменится, если ошибки опустить?

Данные + модель и остатки: модель «не врёт»

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4.2))
t = np.linspace(0, 10, 200)
ax1.scatter(x, y, s=25, label='данные')
ax1.plot(t, model(t, *popt), color='#059669', label='модель')
ax1.set_xlabel('t, с'); ax1.set_ylabel('y, отн. ед.')
ax1.legend(); ax1.grid(alpha=0.3)
ax2.scatter(x, resid, s=25)
ax2.axhline(0, color='gray', linestyle='--')
ax2.set_xlabel('t, с'); ax2.set_ylabel('остатки')
ax2.grid(alpha=0.3)
fig.tight_layout()
fig.savefig('figures/fit_reaction.png', dpi=150)

reaction.csv, 25 точек. Остатки вокруг нуля без тренда и без «воронки» — модель не врёт систематически.

Зачем показывать остатки, если модель «красиво легла» на данные? (глаз обманывает; остатки выявляют систематическую ошибку формы)

Разумный p0: оцениваем по форме данных

# p0 — оценка по форме данных
curve_fit(model, x, y, p0=[1.0, 1.0, 1.0])
# A = 4.94, λ = 0.30, c = 0.53
# RMSE = 0.121   (хорошо)

# p0 = нули — «застрял»
curve_fit(model, x, y, p0=[0.0, 0.0, 0.0])
# A = -19773, λ ≈ 0.0, c = 19777
# RMSE = 0.503   (застрял)

Правило: нелинейная подгонка решает задачу минимизации — она может сойтись в ложном локальном минимуме. Стартовые значения оцениваются по данным: при t=0 значение ≈ A+c ≈ 5.5, «хвост» при t=10 ≈ c ≈ 0.5, скорость спада — λ ≈ 0.3.

Как оценить p0 для модели затухания, не зная истинных параметров? Что подсказывает форма данных?

curve_fit без p0 застрял: A=−19773, c=19777

# без разумного p0 scipy берёт единицы по умолчанию
curve_fit(model, x, y, p0=[0.0, 0.0, 0.0])
RuntimeError: Optimal parameters not found: Number of calls to function
  has reached maxfev = 800.

При неразумных стартовых значениях алгоритм не укладывается в лимит итераций (maxfev=800) — подгонка не сходится.
На других данных/версиях это же выглядит как «мусорные» параметры-гиганты: A=−19773, λ≈0, c=19777, RMSE=0.503 — A и c компенсируют друг друга, RMSE вчетверо хуже.
Диагностика: не сошлось, RMSE большой или параметры абсурдные — проверь p0. Это происходит со всеми; стартовые значения оцениваются по данным, а не угадываются.

Какой симптом неудачи в обоих вариантах — ошибка и «мусорные» параметры? (главное: без разумного p0 результат не заслуживает доверия)

Live-coding 1: полный цикл curve_fit

Задача: воспроизвести цикл подгонки «данные → модель → p0 → fit → ошибки → остатки → график».

Загрузить reaction.csv;
определить model(t, A, lam, c);
оценить p0 по форме данных;
curve_fit;
perr, RMSE;
график «данные + модель» и «остатки»;
savefig.

После блока — «повтори за лектором» 1–2 минуты. Комментируйте вслух: откуда p0, что вернёт pcov, что должно быть на графике остатков.

Что должно быть на графике остатков у «честной» подгонки?

interp1d: линейная и кубическая по 6 узлам

from scipy.interpolate import interp1d

xn = np.array([0, 1, 2, 3, 4, 5])
yn = np.array([0.0, 0.8, 1.5, 1.2, 0.7, 0.3])

f_lin = interp1d(xn, yn, kind='linear')
f_cub = interp1d(xn, yn, kind='cubic')
f_lin(2.5); f_cub(2.5)
1.35  1.44

Правило: интерполяция проходит через узлы и не даёт параметров — в отличие от curve_fit; кубическая — гладкая, линейная — ломаная; в середине интервала они близки (1.35 vs 1.44).

Когда нужен interp1d, а когда curve_fit?

Экстраполяция за диапазон → ValueError

f_lin(6.5)   # вне диапазона узлов
ValueError: A value (6.5) in x_new is above the
interpolation range's maximum value (5).

Правило: по умолчанию interp1d запрещает выход за пределы диапазона — это защита от бессмысленных значений. Чтение ошибки: последняя строка говорит причину — x_new выше максимума.

Почему scipy выбрасывает ошибку вместо того, чтобы молча вернуть значение?

fill_value='extrapolate' → «правдоподобный мусор» −0.3

interp1d(xn, yn, kind='linear',
         fill_value='extrapolate')(6.5)
-0.3   # «мусор»: отрицательное при положительных данных

interp1d в [0,5] проходит через узлы; fill_value='extrapolate' продолжает линию механически — на x=6.5 получается −0.3 при положительных данных.

Все данные положительные — почему на x=6.5 получилось −0.3? (интерполятор не знает физики задачи — он продолжает линию механически)

Когда interp1d, а когда curve_fit

Вопрос	interp1d	curve_fit
Что делает	проходит через узлы	подбирает параметры модели
Выход	функция-интерполятор	параметры ± ошибки, RMSE
За пределами данных	запрещено (мусор)	можно, но с оговоркой границ
Когда использовать	восстановить табличную функцию	есть физическая/эмпирическая модель

Данные сенсора — какой приём? Затухание маятника с известной моделью — какой?

scipy.stats.norm: pdf, cdf, ppf обзорно

from scipy.stats import norm
norm.pdf(0)          # плотность в 0
0.3989422804014327
norm.cdf(1.96)       # P(X ≤ 1.96) ≈ 0.975
0.9750021048517795
norm.ppf(0.975)      # обратная CDF — квантиль
1.959963984540054

Правило: распределения в scipy.stats: pdf — плотность, cdf — функция распределения, ppf — обратная cdf (квантиль). ppf используется для доверительных интервалов: квантиль 0.975 нормального ≈ 1.96.

Как через ppf получить 95% доверительный интервал для стандартной нормальной величины?

ttest_ind: сравнение средних по полу, p=0.271

from scipy.stats import ttest_ind

ages_male = df[df['Sex'] == 'male']['Age'].dropna()
ages_female = df[df['Sex'] == 'female']['Age'].dropna()
ttest_ind(ages_male, ages_female)
TtestResult(statistic=-1.118, pvalue=0.271)

Правило: ttest_ind проверяет гипотезу о равенстве средних двух выборок; p = 0.271 > 0.05 → по нашим данным нет оснований отвергнуть равенство средних. Средние: мужчины 26.7 (n=24), женщины 30.6 (n=14) — различие ~4 года, но статистически незначимо на таких выборках.

Как сформулировать результат корректно: «пол влияет на возраст» или «мы не находим статистически значимой разницы средних»?

pearsonr: Fare–Survived r ≈ −0.22, p ≈ 0.18

from scipy.stats import pearsonr
pearsonr(df['Fare'], df['Survived'])
PearsonRResult(statistic=-0.218, pvalue=0.176)

Правило: pearsonr возвращает (r, p): r — сила и знак линейной связи, p — значимость. r ≈ −0.22, p ≈ 0.18 — связь слабая и статистически незначимая. Это ровно то число, что мы видели в heatmap недели 7.

Можно ли сказать «цена влияет на выживание»? (нельзя — корреляция ≠ причинность; возможен третий фактор, например класс каюты)

Что значит p = 0.271 простыми словами

Правильно: «если бы средние действительно были равны, вероятность получить такое (или более экстремальное) различие средних составила бы ~27%. Это не доказательство отсутствия различия».
Ошибка: «p = вероятность, что гипотеза верна» — так говорить нельзя.
Порог: 0.05 — конвенция, а не закон природы.

Чем различаются две формулировки? Путать это не стыдно — это одна из самых частых ошибок в статистике.

«Взрыв продаж» vs честный рост 22%

Ось Y обрезана (118..128) — рост выглядит как рывок.

Ось от нуля (0..140) — тот же ряд пологий, рост ~22%.

Один и тот же ряд: sales = 100 + 2·месяц + N(0, 2), np.random.default_rng(7).

Что здесь не так? Что скрывает левый график? (реальный рост — 22%, а не «взрыв»)

Что скрывает обрезка: рост 22.3%, а не «рывок»

rng = np.random.default_rng(7)
months = np.arange(1, 13)
sales = 100 + 2 * months + rng.normal(0, 2, 12)
growth = (sales[-1] - sales[0]) / sales[0] * 100
print(f"рост за год: {growth:.1f}%")
рост за год: 22.3%

Реальный рост с 1-го по 12-й месяц — 22.3% (~102 → ~125), то есть ~2% в месяц.
Обрезка оси (118..128) скрывает абсолютный уровень и превращает небольшой тренд в «драматичный».
Обман не в данных — обман в масштабе осей. Маркер разрыва (//) был бы честным способом «зума», но здесь не оправдан.

С каким правилом недели 7 это связано? (бары с нуля)

3D-эффекты, размер маркера, базы процентов

Один и тот же столбец: высота одинаковая, а «масса» и точка чтения верха — разные.

3D/объём у bar: «масса» кажется разной, верх читается по-разному. Честно: плоская 2D-диаграмма.
Размер маркера: площадь круга — квадрат радиуса: двойной радиус = кажущийся рост в 4 раза.
Честные проценты: «вырос на 50%» — от какой базы? всегда уточнять.
Не скрывать выбросы и «нулевые» значения.

Почему «объёмный» столбец «тяжелее» плоского, хотя высота та же?

Принципы честной визуализации

Ось Y с нуля для bar-диаграмм.
Не обрезать ось без маркера разрыва (//).
Не скрывать выбросы.
Не использовать 3D-эффекты «для красоты».
Честные проценты: называть базу.
Подписывать единицы измерения.
Размер маркера пропорционален данным.
Остатки — обязательная часть подгонки модели.

Это требования к коду, а не этикет — на ревью масштаб осей проверяется.

Как проверить каждый пункт одним предложением?

Типичные ошибки недели 8

curve_fit без разумного p0 — застревание в локальном минимуме, параметры «мусорные».
Экстраполяция interp1d за пределы данных — правдоподобные, но бессмысленные значения.
Путаница «подобрать параметры» (curve_fit) и «интерполировать» (interp1d).
Полином высокой степени «для красоты» — переобучение, кривая скачет.
quad с векторной функцией — подынтегральная функция должна быть скалярной.
Остатки не показаны — нельзя оценить, «врёт» ли модель систематически.
Вывод «корреляция → причинность» без оговорок.

Студент подобрал полином 5-й степени, и он «идеально» прошёл через все точки. Почему это плохая модель? (переобучение: на новых точках — огромная ошибка)

Вопросы для проверки понимания

Что возвращает quad и как проверить точность?
Как подобрать параметры экспоненциальной модели к данным с шумом?
Чем интерполяция отличается от экстраполяции и почему последняя опасна?
Какой график «вводит в заблуждение» масштабом и почему?
Чем ttest_ind отличается от pearsonr?

Ответьте письменно за 3 минуты — это мини-самооценка перед семинаром.

Переход: curve_fit — это блок 2 проекта модуля.

Проект модуля: curve_fit — блок 2

Проект модуля 3 (неделя 9) = численный анализ + парсер логов + визуализация в одном ноутбуке.
Этап недели 8 — «curve_fit + честные графики»: модель с обоснованным p0, параметры ± ошибки, остатки, графики данных/модели/остатков.
Требования: без экстраполяции; обоснованный выбор модели; остатки визуализированы; графики подписаны.
Автопроверка: функция fit_model(x, y) с тестами на качество подгонки (RMSE < порога).

Пазл проекта: визуализация (нед. 7) + числа (нед. 8) + текст (нед. 9). Что ваша неделя 8 даёт в финальный ноутбук?

Итоги: 6 навыков и что дальше

Сегодня вы научились:

считать интегралы с оценкой ошибки;
решать ОДУ и сверять с аналитикой;
подбирать параметры через curve_fit с p0/pcov/RMSE;
различать интерполяцию и экстраполяцию;
трактовать p-value и корреляцию честно;
узнавать и чинить манипулирующие графики.

Семинар недели 8: «Интегрирование + статистика» (базовый), «Подбор кривой + интерполяция» (стандартный), «ОДУ + честный график» (продвинутый), «экспонента vs полином» (challenge).
ДЗ недели 8: аппроксимировать экспериментальные данные моделью, обосновать выбор, показать остатки (autograder: fit_model).
One-minute paper: запишите одной фразой главное, что поняли, и один вопрос, что остался.

Откройте семинар-8.md, возьмите reaction.csv и подберите первую модель.