Теоретические выводы с нуля, реализация R и обсуждение байесовской точки зрения
Линейная регрессия может быть установлена и интерпретирована с байесовской точки зрения. Первые части обсуждают теорию и предположения практически с нуля, а последующие части включают реализацию R и замечания. Читатели могут свободно скопировать два блока кода в блокнот R и поэкспериментировать с ним.
Начиная с основ
Напомним, что в линейной регрессии нам задаются целевые значения y, данные X,, и мы используем модель
где y - вектор N * 1, X - матрица N * D, w - вектор D * 1, а ошибка - вектор N * 1. У нас есть N точек данных. Измерение D понимается с точки зрения характеристик, поэтому, если мы используем список x, список x² (и список единиц, соответствующих w_0), мы говорим D = 3. Если вам не нравится матричная форма, подумайте о ней как о сокращенной форме следующего, где все является масштабирующим устройством, а не вектором или матрицей:
В классической линейной регрессии предполагается, что член ошибки имеет нормальное распределение, поэтому сразу следует, что y нормально распределено со средним значением Xw, а дисперсия любой дисперсии - ошибка термин имеет (обозначается σ² или диагональной матрицей с элементами σ²). Нормальное предположение в большинстве случаев оказывается правильным, и эту нормальную модель мы также используем в байесовской регрессии.
Вывод параметров
Теперь мы сталкиваемся с двумя проблемами: вывод w и предсказание y для любого нового X. Используя хорошо известное правило Байеса и приведенные выше предположения, мы всего лишь на шаг впереди не только к решению этих двух проблем, но и к получению полное распределение вероятностей y для любого нового X. Вот правило Байеса с использованием наших обозначений, которое выражает апостериорное распределение параметра w для заданных данных:
π и f - функции плотности вероятности. Поскольку результат является функцией от w, мы можем игнорировать знаменатель, зная, что числитель пропорционален левой части на константу. Мы знаем из предположений, что функция правдоподобия f (y | w, x) следует нормальному распределению. Другой термин - это априорное распределение w, и это отражает, как следует из названия, априорное знание параметров.
Предварительное распределение. Определение предшествующего распределения - интересная часть байесовского рабочего процесса. Для удобства мы положим w ~ N (m_o, S_o), а гиперпараметры m и S теперь отражают предварительное знание w. Если вы плохо разбираетесь в w или находите какое-либо присвоение m и S слишком субъективным, поправкой можно считать «неинформативные» априорные значения. В этом случае мы устанавливаем m равным 0 и, что более важно, устанавливаем S как диагональную матрицу с очень большими значениями. Мы говорим, что w имеет очень высокую дисперсию, и поэтому мы мало знаем, каким будет w.
После определения всех этих функций вероятности несколько строк простых алгебраических манипуляций (на самом деле довольно много строк) дадут апостериорную оценку после наблюдения N точек данных:
Это похоже на набор символов, но все они уже определены, и вы можете вычислить это распределение, как только этот теоретический результат будет реализован в коде. (N (m, S) означает нормальное распределение со средним m и ковариационной матрицей S.)
Прогнозирующее распространение
Полный байесовский подход означает не только получение единственного прогноза (обозначьте новую пару данных y_o, x_o), но также получение распределения этой новой точки.
Мы сделали обратное маргинализации от соединения, чтобы получить маржинальное распределение в первой строке, и использование правила Байеса внутри интеграла во второй строке, где мы также удалили ненужные зависимости. Обратите внимание, что нам известны две последние функции вероятности. Результат полного прогнозного распределения:
Реализация в R
Реализация в R довольно удобна. Опираясь на приведенные выше теоретические результаты, мы просто вводим матричные умножения в наш код и получаем результаты как прогнозов, так и прогнозных распределений. Чтобы проиллюстрировать это примером, мы используем игрушечную задачу: X равно от -1 до 1, равномерно, а y строится как следующие дополнения синусоидальных кривых с нормальным шумом (см. График ниже для иллюстрации y).
Следующий код получает эти данные.
library(ggplot2) # — — — — — Get data — — — — — — — — — — — — — — — — — — — — —+ X <- (-30:30)/30 N <- length(X) D <- 10 var <- 0.15*0.15 e <- rnorm(N,0,var^0.5) EY <- sin(2*pi*X)*(X<=0) + 0.5*sin(4*pi*X)*(X>0) Y <- sin(2*pi*X)*(X<=0) + 0.5*sin(4*pi*X)*(X>0) + e data <- data.frame(X,Y) g1 <- ggplot(data=data) + geom_point(mapping=aes(x=X,y=Y))
Следующий код (в разделе «Вывод») реализует вышеуказанные теоретические результаты. Мы также расширяем возможности x (обозначается в коде как phi_X в разделе «Создание базовых функций»). Так же, как мы расширили бы x до x² и т. Д., Теперь мы расширяем его до 9 радиальных базисных функций, каждая из которых выглядит следующим образом. Обратите внимание, что хотя они выглядят как нормальная плотность, они не интерпретируются как вероятности.
Одним из преимуществ радиальных базисных функций является то, что радиальные базисные функции могут соответствовать множеству кривых, включая полиномиальные и синусоидальные.
# — — — — — Construct basis functions — — — — — — — — — — — —+
phi_X <- matrix(0, nrow=N, ncol=D)
phi_X[,1] <- X
mu <- seq(min(X),max(X),length.out=D+1)
mu <- mu[c(-1,-length(mu))]
for(i in 2:D){
phi_X[,i] <- exp(-(X-mu[i-1])^2/(2*var))
}
# — — — — — Inference — — — — — — — — — — — — — — — — — — — —+
# Commented out is general prior
# m0 <- matrix(0,D,1)
# S0 <- diag(x=1000,D,D)
# SN <- inv(inv(S0)+t(phi_X)%*%phi_X/var)
# mN <- SN%*%(inv(S0)%*%m0 + t(phi_X)%*%Y/var)
# Y_hat <- t(mN) %*% t(phi_X)
# We use non-informative prior for now
m0 <- matrix(0,D,1)
SN <- solve(t(phi_X)%*%phi_X/var)
mN <- SN%*%t(phi_X)%*%Y/var
Y_hat <- t(mN) %*% t(phi_X)
var_hat <- array(0, N)
for(i in 1:N){
var_hat[i] <- var + phi_X[i,]%*%SN%*%phi_X[i,]
}
g_bayes <- g1 +
geom_line(mapping=aes(x=X,y=Y_hat[1,]),color=’#0000FF’)
g_bayes_full <- g_bayes + geom_ribbon(mapping=aes(x=X,y=Y_hat[1,],
ymin=Y_hat[1,]-1.96*var_hat^0.5,
ymax=Y_hat[1,]+1.96*var_hat^0.5, alpha=0.1),
fill=’#9999FF’)
Одна деталь, которую следует отметить в этих вычислениях, заключается в том, что мы используем неинформативный априор. Закомментированный раздел - это в точности теоретические результаты, приведенные выше, в то время как для неинформативных предварительных результатов мы используем ковариационную матрицу с диагональными элементами, приближающимися к бесконечности, поэтому обратная величина в этом коде напрямую рассматривается как 0. Если вы хотите использовать этот код, убедитесь, что вы установили пакет ggplot2 для построения графика.
Следующая иллюстрация предназначена для представления полного прогнозного распределения и дает представление о том, насколько хорошо данные соответствуют.
Замечания
Результат множественной линейной регрессии такой же, как и в случае байесовской регрессии с использованием неправильного априорного значения с бесконечной ковариационной матрицей. Как правило, хорошей практикой является получение некоторых эмпирических знаний о параметрах и использование информативного априорного значения. Затем байесовская регрессия может количественно оценить и показать, как различные предшествующие знания влияют на прогнозы. В любом случае байесовское представление может легко интерпретировать диапазон прогнозов y как вероятность, отличную от доверительного интервала, вычисленного с помощью классической линейной регрессии.
Подбор данных с этой точки зрения также упрощает «обучение на ходу». Скажем, я сначала наблюдал 10000 точек данных и вычислил апостериорную величину параметра w. После этого мне каким-то образом удалось получить еще 1000 точек данных, и вместо того, чтобы снова запускать всю регрессию, я могу использовать ранее вычисленные апостериорные данные в качестве априорных для этих 1000 точек. Этот последовательный процесс дает тот же результат, что и повторное использование всех данных. Мне нравится эта идея, поскольку она интуитивно понятна, так как усвоенное мнение пропорционально ранее усвоенным мнениям плюс новым наблюдениям, а обучение продолжается. В анекдоте говорится, что байесовец, который мечтает о лошади и наблюдает за ослом, назовет ее мулом. Но если он сделает больше наблюдений за ним, в конце концов он скажет, что это действительно осел.
