Перейдем к оценке байесовской регрессии пик-плато.
Для этого мы воспользуемся встроенным набором
данных по машинам: регрессия, spike and slab regression.
Соответственно, мы в набор данных h поместим данные по машинам.
Ну давайте для того, чтобы работать не с милями и футами,
а с километрам и метрами, мы сначала его изменим.
То есть мы сделаем mutate, изменим набор данных cars, и формулы,
по которым мы будем менять скорость в километрах в час — это 1,61 умножить
на скорость в милях в час исходное, которое была в массиве данных и,
соответственно, длина тормозного пути, чтобы перевести ее
из футов в метры надо помножить на 0,3 исходную длину в наборе данных cars.
Итого, вот у нас есть привычные переменные,
давайте мы в этот набор данных h, во-первых, добавим переменную, которую
назовем мусор — просто бессмысленная переменная чтобы проэспериментировать и
посмотреть как байесовская регрессия относится к бессмысленной переменной.
Значит, здесь мы просто создадим, сгенерим
нормальные стандартные случайные величины.
Это будет искусственная мусорная переменная.
Также давайте в набор данных мы добавим переменную квадрат,
мы изменим набор данных h и туда добавим квадрат скорости,
ну который, естественно, будет равен скорости в квадрате.
И теперь мы построим обычную модель,
модель метода наименьших квадратов по набору h.
Мы построим регрессию длины тормозного
пути на скорость и мусорную переменную, которая просто так вносит шум.
Здесь плюс конечно надо поставить, плюсик.
Вот и давайте посмотрим на результаты оценивания классической линейной модели.
То есть мы видим, что здесь у нас коэффициент, конечно, есть и при скорости
и при переменной мусор, но переменная мусор она у нас статистически незначима.
А теперь мы по тем же самым данных оценим регрессию пик-плато то же самой модели.
И на этот раз мы сможем задать прямой вопрос: а чему же равна вероятность того,
что коэффициент при мусорной переменной на самом деле в точности ноль.
model_ spike and
slab regression равняется.
Функция называется spike and slab.
Ну здесь надо отметить, что на самом деле функция spike and slab оценивает,
если быть честным, гораздо более сложные модели, и это некоторое упрощение
взгляда того, что там только возможен коэффициент ноль или не ноль,
там самом деле немножко более тонкий вариант, но примерно такой.
Соответственно, мы возьмем данные из набора данных h,
формула по которой мы оцениваем будет точно такая же,
то есть мы изучаем зависимость длины тормозного пути от истинной скорости
машины и некоторой искусственной мусорной переменной,
и давайте зададим размер выборки.
По умолчанию, по-моему, алгоритм генерит 500 наблюдений из
апостериорного распределения, но тут чем больше точность,
чем больше наблюдений из апостериорного распределения, тем выше точность.
И единственное, чем вы ограничены — это временем вычислений.
То есть чем больше выборка из апостериорного вычисления...
распределения, тем точнее результат.
Ну разве что вам нужно получить результат очень быстрый, а увеличение выборки
в 2 раза приводит к увеличению времени там до двух дней или там до трех.
Но в данном случае такого не произойдет, поэтому мы просто возьмем, скажем,
4000 наблюдений.
И, соответственно, запустим оцениваться нашу модель.
Ну она занимает долгое время, потому что алгоритм MCMC он не просто минимизирует,
а там генерятся случайные величины, и на каждом шагу много операций.
Можем посмотреть, соответственно, результаты оценивания.
Поскольку пакет по spike and slab regression еще не очень доработанный,
то тут немножко нестандартный синтаксис: не summary, а print model_ss.
И, соответственно, мы видим,
что здесь автоматом тоже отобралась переменная speed,
но мы можем более точно заглянуть что делал пакет,
можем посмотреть model_ss$ summary.
И вот здесь вот можем увидеть все-таки оценки
коэффициентов, которые были за кадром.
Вот bma.scale — это оценки коэффициентов в оригинальном масштабе.
Дело в том, что пакет на самом деле при осуществлении своих операций он
сначала отмасштабировал все переменные.
Вот. Ну и здесь видно,
что значит с увеличением скорости на единичку длина тормозного
пути растет на 0,69, и, соответственно, мусорная переменная,
с увеличением ее на единичку, длина тормозного пути растет на 0,20.
Но, тем не менее, алгоритм автоматом убрал мусорную переменную.
Ну и можно все-таки посмотреть на вероятности.
Давайте мы посмотрим какие коэффициенты включались в
каждую их 4000 случайно оцененных,
выборку из 4000 апостериорных значений гамма.
То есть гамма — это индикатор включалась или не включалась переменная.
Соответственно, давайте мы введем такой массив данных included_regressors,
и мы возьмем из
модели ss список,
который называется model.
И применим функцию melt,
которая превратит это в привычный набор данных, в привычный dataframe.
Ну и давайте чтобы понять, что такое included
regressors просто посмотрим, что это такое мы страшное создали.
Это такой здоровый список,
в котором для каждого значения выборки из апостериорного распределения
указано какие регрессоры туда включались.
Ну поскольку модель случайным способом оценивается,
то у вас при симуляции могут получаться другие результаты.
Ну, например, я про свои расскажу.
Вот ежели я возьму 3031-е
значение выборки из апостериорного распределения,
то там был и первый регрессор, и второй.
Если я возьму, скажем, 3020-е значение выборки из апостериорного распределения,
то там включался только первый регрессор, то там включался только первый регрессор.
Ну и соответственно, поскольку я знаю, что всего у меня размер 4000, то я могу,
ну давайте я посмотрю на начало набора данных included_regressors вот,
соответственно, здесь есть value и l1.
l1 — это номер наблюдений в выборке,
а value — это какие регрессоры туда включались.
Вот видно, что в 1-ю,
2-ю и в 3-ю и в 4-ю значения выборки входил только первый регрессор.
И я могу, например, посчитать сумму,
included_reressors$value равняется,
равняется единичке, и поделить на 4000 — общее количество элементов в выборке,
соответственно, у меня получается, что регрессор скорость включался в
каждую из значений апостериорного распределения.
А если я посмотрю, какова вероятность того,
что переменная мусорная входит в модель, у нее номер 2 и,
соответственно, здесь я получу вероятность 0,35.
Соответственно, согласно байесовскому подходу, у меня получается,
что с вероятностью 1 переменная скорость
влияет на, соответственно, длину тормозного пути.
Ну и поскольку у нас маленькая выборка, в нашей выборке всего,
наша выборка называется h nrow(h),
в нашей выборке всего 50 наблюдений, соответственно, по нашей выборке,
когда мы добавили искусственную переменную, из-за которой был только шум,
ну по имеющимся данным, при имеющемся априорном распределении, вероятность того,
что мусорная переменная влияет на длину тормозного пути равна 0,35.
И на этом примере байесовской
регрессии пик-плато мы закончим наши компьютерные занятия в этот раз.