Байесовская линейная регрессия

Обычная регрессия выдаёт одну линию. Байесовская оставляет целое облако линий, взвешенных по тому, насколько каждая согласуется и с данными, и с тем, что мы знали до них. Там, где измерений мало, облако расходится веером — и честно расширяет прогнозный интервал ровно настолько, насколько мы не знаем ответа.

Шесть вечеров и целый веер прямых

Возьмём реальные данные проката велосипедов: вечерний час пик, 17:00, только рабочие дни 2011 года — всего 250 вечеров. Признак один: нормированная температура temp. Ответ — число поездок за час. По всем 250 вечерам метод наименьших квадратов из урока 49 даёт наклон 599 поездок на полный размах температуры и средний уровень 395 поездок; типичный разброс вокруг прямой — 111 поездок.

Теперь представим, что мы — лаборатория, у которой измерений почти нет. Возьмём случайные шесть вечеров из этих же 250. МНК по шести точкам выдаёт наклон 799. Ошибка относительно «полного» наклона — 200 поездок, треть величины. И главное: МНК не сообщает об этом ни слова. Он возвращает число, как если бы оно было измерено, а не угадано по шести шумным точкам.

Байесовский подход отвечает иначе. Он не выбирает прямую — он приписывает каждой прямой вес и держит их все. Через область шести точек проходит плотный пучок; левее и правее пучок расходится. Наклон получается 542 с собственным стандартным отклонением 170 — и, что приятно, эта оценка ближе к «полному» наклону 599, чем МНК: промах 57 вместо 200.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Что мы считаем неизвестным

В уроке 49 вектор весов $w$ был неизвестной константой, которую мы оцениваем. Здесь мы делаем шаг из урока 47 и объявляем $w$ случайным вектором: у него есть распределение до данных и распределение после. Модель наблюдений остаётся прежней:

$$y=Xw+\varepsilon,\qquad \varepsilon\sim N(0,\sigma^2I),$$

то есть

$$p(y\mid X,w)=\prod_{i=1}^{n} \frac{1}{\sqrt{2\pi}\,\sigma} \exp\left[-\frac{(y_i-x_i^\top w)^2}{2\sigma^2}\right].$$

К ней добавляется нормальный prior

$$w\sim N(m_0,S_0),$$

а правило Байеса из урока 42 в форме плотностей даёт

$$p(w\mid D)\;\propto\;p(y\mid X,w)\,p(w).$$

Prior — это не «вера в удачу», а записанное предметное знание. Для велопроката мы взяли $m_0=0$ и $S_0=\tau^2I$ с $\tau=300$: заранее считаем, что переход от самого холодного вечера к самому тёплому вряд ли меняет спрос больше чем на несколько сотен поездок, а знак эффекта не навязываем. Это слабое, но не пустое утверждение — и именно оно удержало наклон от прыжка на 799.

Две нормальные формы умножаются в третью

Логарифм правдоподобия по $w$ — квадратичная форма:

$$\ln p(y\mid X,w)= -\frac{1}{2\sigma^2}(y-Xw)^\top(y-Xw)+\text{const},$$

логарифм prior — тоже:

$$\ln p(w)=-\frac12(w-m_0)^\top S_0^{-1}(w-m_0)+\text{const}.$$

Сумма двух квадратичных форм — снова квадратичная форма, а экспонента от квадратичной формы — снова нормальная плотность. Значит, posterior нормален, и достаточно собрать коэффициенты. При $m_0=0$ показатель равен

$$-\frac12\,w^\top\!\left(S_0^{-1}+\frac{X^\top X}{\sigma^2}\right)\!w +w^\top\frac{X^\top y}{\sigma^2}+\text{const},$$

откуда, сравнивая с общим видом $-\tfrac12 w^\top S_n^{-1}w+w^\top S_n^{-1}m_n$,

$$S_n^{-1}=S_0^{-1}+\frac{1}{\sigma^2}X^\top X,$$ $$m_n=S_n\left(S_0^{-1}m_0+\frac{1}{\sigma^2}X^\top y\right),$$ $$w\mid D\sim N(m_n,S_n).$$

Складываются не ковариации, а точности — матрицы, обратные ковариациям. Prior приносит точность $S_0^{-1}$, данные — точность $X^\top X/\sigma^2$. Это дословно многомерная версия сопряжённого обновления из урока 48, где складывались счётчики бета-распределения.

Обратите внимание на асимметрию: $S_n$ не зависит от $y$. Ширина posterior определяется тем, где мы измеряли и насколько шумен прибор, но не тем, что получилось. Это свойство нормальной модели, и оно же делает возможным планирование эксперимента заранее — к нему мы вернёмся ниже.

Одно измерение на бумаге

Проследим формулы на числах, которые можно проверить в уме. Пусть модель одномерна: $y=xw+\varepsilon$, шум $\varepsilon\sim N(0,4)$, prior $w\sim N(0,1)$. Получено одно измерение: $x=2$, $y=6$. Точность posterior:

$$S_n^{-1}=\frac{1}{1}+\frac{x^2}{\sigma^2}=1+\frac{4}{4}=2, \qquad S_n=0{,}5 .$$

Среднее:

$$m_n=S_n\cdot\frac{xy}{\sigma^2}=0{,}5\cdot\frac{2\cdot 6}{4}=1{,}5 .$$

МНК по одной точке дал бы $w=y/x=3$. Prior и данные принесли одинаковую точность (по единице каждый), поэтому ответ встал ровно посередине между нулём prior и тройкой данных. Никакой мистики: посередине — потому что голоса равны.

Для нового входа $x_*=3$ средний прогноз равен $4{,}5$. Дисперсия самой линии составляет $x_*^2S_n=9\cdot0{,}5=4{,}5$, а будущее наблюдение добавляет ещё $\sigma^2=4$:

$$\operatorname{Var}(y_*\mid x_*,D)=4{,}5+4=8{,}5,\qquad \sqrt{8{,}5}\approx 2{,}92 .$$

Сообщить одно только $4{,}5$ значит потерять обе причины разброса сразу.

Ridge — это posterior с нормальным prior

Возьмём $m_0=0$, $S_0=\tau^2I$. Тогда

$$m_n=\left(X^\top X+\frac{\sigma^2}{\tau^2}I\right)^{-1}X^\top y,$$

то есть в точности решение гребневой регрессии из урока 51 с параметром

$$\lambda=\frac{\sigma^2}{\tau^2}.$$

Это не аналогия, а тождество. На реальном наборе diabetes (442 пациента, 10 стандартизованных признаков, оценка шума $\widehat\sigma=54{,}1$) мы прогнали шестьдесят значений $\tau$ и сравнили posterior mean с решением ridge: обе формулы совпали до машинной точности, максимальное расхождение по всем коэффициентам и всем $\tau$ оказалось меньше $10^{-8}$.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Разница между ridge и байесовской моделью не в точке ответа, а в том, что байесовская модель отдаёт ещё и $S_n$. Ridge говорит «наклон 542». Posterior говорит «наклон 542 плюс-минус 170» — и это совсем другое сообщение руководителю проекта.

Гаусс вывел метод наименьших квадратов именно как поиск наиболее вероятного значения при нормальных ошибках — то есть, на нашем языке, как поиск моды posterior при плоском prior. Байесовская регрессия просто отказывается выбрасывать всё, кроме моды.

Эллипс незнания: сумма известна, разность — нет

В наборе велопроката есть два почти одинаковых признака: temp (температура) и atemp (ощущаемая температура). Их корреляция на наших 250 вечерах равна $0{,}991$. Что произойдёт, если подать в модель оба?

Правдоподобие вытянется в узкий гребень: данные хорошо определяют сумму коэффициентов (потому что признаки почти совпадают и работает именно сумма) и почти не определяют их разность. Prior круглый, он ограничивает общую норму вектора. Произведение даёт эллипс покороче гребня, но всё ещё вытянутый.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Числа таковы: стандартное отклонение posterior вдоль направления суммы $(w_1+w_2)/\sqrt2$ равно 22 поездкам, а вдоль направления разности $(w_1-w_2)/\sqrt2$ — 200 поездок, в 9,1 раза больше. Сумма коэффициентов получилась $129$ и определена уверенно; сами коэффициенты по отдельности — почти нет.

Предсказательное распределение: два слоя незнания

Для нового объекта $x_*$ прогноз при фиксированных весах равен $x_*^\top w$. Усредним по posterior. Линейная функция нормального вектора нормальна, поэтому

$$x_*^\top w\mid D\sim N\!\left(x_*^\top m_n,\;x_*^\top S_nx_*\right),$$

а будущее наблюдение добавляет собственный шум:

$$y_*\mid x_*,D\sim N\!\left(x_*^\top m_n,\;\sigma^2+x_*^\top S_nx_*\right).$$

Формально это интеграл

$$p(y_*\mid x_*,D)=\int N(y_*\mid x_*^\top w,\sigma^2)\,N(w\mid m_n,S_n)\,dw,$$

и он берётся в явном виде именно потому, что обе части нормальны. Дисперсия распадается на два слагаемых с разной судьбой:

$$\operatorname{Var}(y_*\mid x_*,D)= \underbrace{\sigma^2}_{\text{шум мира}}+ \underbrace{x_*^\top S_nx_*}_{\text{незнание весов}} .$$

Второе слагаемое стремится к нулю с ростом числа наблюдений — незнание поправимо. Первое не исчезает никогда: даже зная веса точно, мы не знаем, сколько именно человек завтра решит поехать. Это ровно то различие между интервалом для среднего и интервалом для нового наблюдения, которое мы разбирали в уроке 46.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Почему веер раскрывается на краю

Запишем предсказательную дисперсию для одномерной модели с центрированным признаком, $\phi(x)=(1,x-\bar x)$:

$$x^\top S_nx=S_{00}+2(x-\bar x)S_{01}+(x-\bar x)^2S_{11}.$$

Это парабола по $x$ с минимумом вблизи центра данных. Геометрически всё просто: прямые из posterior поворачиваются вокруг области наблюдений, как ножницы, и чем дальше от точки опоры, тем больше расхождение концов. На наших шести вечерах стандартное отклонение среднего равно 45 поездкам при temp $=0{,}5$ и 98 поездкам при temp $=0$ — веер раскрылся в 2,2 раза.

Лаборатория веера

Внутри — те же данные: вечера велопроката из той же выборки 2011 года, только шестёрка взята другая, поэтому и числа будут другие. Начните с шести точек и включите режим «полосы и линии»: сорок прямых из posterior наглядно показывают, что именно скрыто за полосой. Потом растяните ползунок ширины prior: при $\tau=25$ наклон почти прижат к нулю (модель отказывается верить, что температура вообще влияет), при $\tau=800$ синяя линия почти ложится на красную штриховую линию МНК — почти, а не точно: prior конечной ширины никогда не отпускает оценку до конца. Увеличьте шум $\sigma$: золотая полоса раздувается всюду, а синяя — только вслед за ростом $x^\top S_nx$. Наконец переключитесь на 12 вечеров и заметьте, где веер схлопнулся, а где остался: он остаётся там, где нет точек. Щёлкните по пустому полю слева — вы добавили измерение на холодном краю, и левый край веера сжался сильнее, чем от точки, добавленной в центре.

Когда σ неизвестна: хвосты Стьюдента

До сих пор $\sigma^2$ считалась известной. На практике её оценивают по тем же данным. Сопряжённый выбор — нормально-обратно-гамма prior:

$$\sigma^2\sim\mathrm{InvGamma}(a_0,b_0),\qquad w\mid\sigma^2\sim N(m_0,\sigma^2V_0).$$

Posterior остаётся в том же семействе, а предсказательное распределение, проинтегрированное по $\sigma^2$, оказывается не нормальным, а стьюдентовским:

$$y_*\mid x_*,D\;\sim\; t_{2a_n}\!\left(x_*^\top m_n,\; \frac{b_n}{a_n}\left(1+x_*^\top V_nx_*\right)\right).$$

Смысл прост: сомнение в масштабе шума утяжеляет хвосты. Квантиль 97,5% для нормального распределения равен 1,96, для $t$ с тремя степенями свободы — 3,18. Подставить одну оценку $\widehat\sigma$ и взять 1,96 — значит сделать вид, что масштаб известен точно.

Проверка обещания: калибровка интервала

Интервал, который называет себя 95-процентным, обязан накрывать истину в 95 случаях из ста. Это проверяемое утверждение — и его надо проверять. Здесь мы переходим к явно модельному эксперименту с фиксированным seed: настоящих данных с известной истиной не бывает, а нам нужна именно известная истина.

Схема: истинная прямая $y=1+2x$, шум $\sigma=2$, всего $n=5$ наблюдений, 4000 повторов. В каждом повторе строим 95-процентный предсказательный интервал для нового объекта тремя способами и смотрим, попал ли он.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Вывод общий и выходит далеко за рамки регрессии: если неизвестная величина влияет на прогноз, её неопределённость нужно либо проинтегрировать, либо обоснованно показать, что ею можно пренебречь. Подстановка точечной оценки — это молчаливое утверждение «я знаю это точно», и цена такого утверждения измеряется в недопокрытии.

Куда поставить седьмое измерение

Раз $S_n$ не зависит от $y$, можно спросить: если разрешено одно новое измерение, где его сделать? Добавление точки $x_{n+1}$ меняет точность на матрицу ранга один:

$$S_{n+1}^{-1}=S_n^{-1}+\frac{1}{\sigma^2}\,\phi\phi^\top, \qquad \phi=\phi(x_{n+1}),$$

а по формуле Шермана–Моррисона

$$S_{n+1}=S_n-\frac{S_n\phi\phi^\top S_n}{\sigma^2+\phi^\top S_n\phi}.$$

Из неё видно: сокращение неопределённости тем больше, чем больше $\phi^\top S_n\phi$ — то есть чем меньше модель уверена в этой точке. Критерии выбора называют буквами: D-оптимальность минимизирует $\det S_{n+1}$ (объём эллипсоида незнания), A-оптимальность — след, G-оптимальность — максимум предсказательной дисперсии по области.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

И тут начинается взрослая часть. Точка с максимальной математической неопределённостью может быть недоступна, дорога или опасна: temp $=0$ — это мороз, в который поездок почти не бывает, а при калибровке химического датчика предельная концентрация может быть попросту ядовитой. Кроме того, D-оптимальный план ставит точки только на края — и потому не способен обнаружить кривизну: по двум кучкам на концах прямая и парабола неразличимы. Оптимальность всегда оптимальна относительно принятой модели.

Валерий Фёдоров и теория оптимального эксперимента

Идея «модель сама подсказывает, где мерить дальше» превратилась в математическую дисциплину в 1960–1970-е годы, и большой вклад в неё внесла московская школа. Валерий Вадимович Фёдоров, работавший тогда на кафедре математической статистики МГУ, в 1971 году выпустил монографию «Теория оптимального эксперимента» — книгу, которая уже через год вышла по-английски и стала одним из самых цитируемых текстов по планированию эксперимента в мире.

Фёдоров придал теории вычислительную форму: он разработал последовательные алгоритмы построения оптимальных планов, в которых точки добавляются по одной — каждая туда, где текущий план даёт наибольшую предсказательную дисперсию. Это буквально та процедура, которую мы только что проделали на правой панели рисунка 52.6, и её сходимость к оптимальному непрерывному плану опирается на теорему эквивалентности Кифера–Вольфовица: план D-оптимален тогда и только тогда, когда максимум $\phi^\top M^{-1}\phi$ по области равен числу параметров.

Для нас важно следствие: активное обучение, о котором сегодня говорят в связи с разметкой данных для нейросетей, — прямой потомок этой теории. Модель не только отвечает на вопросы, но и предлагает следующий вопрос миру, и предлагает его, глядя на собственную ковариацию.

От весов к функциям: шаг к гауссовскому процессу

Prior на веса индуцирует распределение на самих функциях. Если $w\sim N(0,S_0)$, то для любых двух точек

$$\operatorname{Cov}\big(f(x),f(x')\big)= \mathbb E\big[x^\top w\,w^\top x'\big]=x^\top S_0x' .$$

Заменив исходную координату богатым набором базисных функций из урока 50, получаем

$$k(x,x')=\phi(x)^\top S_0\,\phi(x'),$$

и предсказательные формулы можно переписать через одну лишь функцию $k$. Гауссовский процесс делает следующий шаг: задаёт $k$ напрямую, не перечисляя базис — иногда бесконечный. Модель тогда описывается не «какие признаки», а «какие точки считать соседними»:

$$k_{\text{RBF}}(x,x')=\exp\!\left[-\frac{(x-x')^2}{2\ell^2}\right], \qquad k_{\text{per}}(x,x')=\exp\!\left[-\frac{2\sin^2\big(\pi|x-x'|/p\big)}{\ell^2}\right].$$

Периодическое ядро связывает одинаковые часы разных суток, гладкое радиальное — близкие координаты, линейное — похожие направления. И снова вывод тот же, что в уроке 50: неопределённость зависит не от мира, а от выбранного представления.

Что писать в отчёте

Полезный отчёт по байесовской регрессии содержит четыре вещи. Первое — posterior mean коэффициентов с их совместной неопределённостью, а не таблицу голых чисел: если два признака коррелируют, отдельные стандартные отклонения обманывают, и лучше показать ковариацию или несколько прямых из posterior. Второе — предсказательные интервалы, явно разделённые на интервал для среднего и интервал для нового объекта. Третье — описание prior на предметном языке («эффект температуры вряд ли превышает несколько сотен поездок») вместе с анализом чувствительности: как изменится вывод при $\tau$ вдвое больше и вдвое меньше. Четвёртое — проверка воспроизведения структуры данных: сгенерируйте выборки из предсказательного распределения и сравните их гистограмму с реальной.

Чего интервал не знает

Байесовская линейная регрессия честна ровно в тех пределах, в которых честны её допущения. Она предполагает линейность по весам, нормальный шум с постоянной дисперсией, независимость наблюдений и правильный набор признаков. Каждое из этих допущений проверяемо: кривизна видна на графике остатков, растущий разброс требует $\sigma(x)$ или логарифма ответа, зависимость наблюдений во времени ломает формулу $X^\top X/\sigma^2$, а забытый признак не выдаст себя вовсе — его эффект просто впитается в шум и в веса соседей.

И всё же главное достижение этой главы — смена формы ответа. Мы вышли из урока с распределением вместо числа. Наклон 542 плюс-минус 170; прогноз на тёплый вечер 565 поездок с интервалом от 297 до 833; полоса, которая расширяется там, где мы не смотрели. Такой ответ труднее вставить в презентацию и невозможно подать как окончательный, но именно он позволяет спросить: достаточно ли мы знаем, чтобы действовать? В уроке 53 мы посмотрим, что делать, когда сопряжённых формул нет и posterior приходится добывать численно.

Задачи