Оценить апостериорное прогнозирующее распределение в байесовской линейной регрессии

10

Я запутался в том, как оценивать апостериорное предиктивное распределение для байесовской линейной регрессии, за пределами основного случая, описанного здесь на странице 3 и скопированного ниже.

p(y~y)=p(y~β,σ2)p(β,σ2y)

Основной случай - это модель линейной регрессии:

y=Xβ+ϵ,yN(Xβ,σ2)

Если мы используем либо униформу, предшествующую , с масштабом-Inv предшествующим , ИЛИ нормальным-обратным гамма-предшествующим (см. Здесь ), то апостериорное предиктивное распределение является аналитическим и является учеником t. χ 2 σ 2βχ2σ2

Как насчет этой модели?

y=Xβ+ϵ,yN(Xβ,Σ)

Когда , но известен, апостериорное предиктивное распределение является многомерным гауссовским. Обычно вы не знаете , но должны оценить это. Может быть, вы скажете его диагональ и каким-то образом сделаете диагональ функцией ковариат. Это обсуждается в главе о линейной регрессии байесовского анализа данных Гельмана .Σ ΣyN(Xβ,Σ)ΣΣ

Существует ли аналитическая форма для апостериорного прогнозирующего распределения в этом случае? Могу ли я просто включить мою оценку этого в многомерный студент т? Если вы оцениваете более одной дисперсии, является ли распределение по-прежнему многовариантным?

Я спрашиваю , потому что у меня есть некоторые уже на руках. Я хочу знать, было ли это более вероятно предсказано, например, линейной регрессией A, линейной регрессией B y~

bill_e
источник
1
Если у вас есть задние выборки из заднего распределения, вы можете оценить прогнозирующее распределение с помощью аппроксимации Монте-Карло
niandra82
Ах, спасибо, я всегда мог это сделать. Нет ли аналитической формулы в этом случае?
bill_e
Ссылки битые, кстати. Было бы здорово, если бы вы включили ссылки по-другому.
Maxim.K

Ответы:

6

Если вы предполагаете униформу до , то для апостериор будет с Чтобы найти прогнозирующее распределение, нам нужно больше информации. Если и условно не зависит от заданного , то Но обычно для этих типов моделей и не являются условно независимыми, вместо этого мы обычно имеем ββ

β|yN(β^,Vβ).
β^=[XΣ1X]XyandVβ=[XΣ1X]1.
y~N(X~β,Σ~)yβ
y~|yN(X~β^,Σ~+Vβ).
yy~~У| у~N(~Хβ+Σ21Σ-1(у-Хβ),~Σ-Σ21Σ-1Σ12). Σ,Σ12,˜Σ
(yy~)N([XβX~β],[ΣΣ12Σ21Σ~]).
Если это так, то Это предполагает, что и все известны. Как вы указываете, как правило, они неизвестны и должны быть оценены. Для общих моделей с такой структурой, например, временных рядов и пространственных моделей, обычно не существует закрытой формы для прогнозирующего распределения.
y~|yN(X~β^+Σ21Σ1(yXβ^),Σ~Σ21Σ1Σ12).
Σ,Σ12,Σ~
jaradniemi
источник
2

При неинформативном или многовариантном априоре Normal-Wishart у вас есть аналитическая форма в виде многовариантного распределения Стьюдента для классической многовариантной множественной регрессии. Я думаю, что события в этом документе связаны с вашим вопросом (вам может понравиться Приложение A :-)). Обычно я сравнивал результаты с последующим прогнозным распределением, полученным с помощью WinBUGS, и аналитической формой: они в точности эквивалентны. Проблема становится сложной только тогда, когда у вас есть дополнительные случайные эффекты в моделях со смешанными эффектами, особенно в несбалансированном дизайне.

В общем случае с классическими регрессиями y и ỹ условно независимы (остатки iid)! Конечно, если это не так, то предлагаемое решение здесь не является правильным.

В R (здесь решение для единообразных априоров), если вы сделали модель lm (названную «модель») одного из откликов в вашей модели и назвали ее «моделью», здесь описано, как получить многомерное прогнозирующее распределение

library(mvtnorm)
Y = as.matrix(datas[,c("resp1","resp2","resp3")])
X =  model.matrix(delete.response(terms(model)), 
           data, model$contrasts)
XprimeX  = t(X) %*% X
XprimeXinv = solve(xprimex)
hatB =  xprimexinv %*% t(X) %*% Y
A = t(Y - X%*%hatB)%*% (Y-X%*%hatB)
F = ncol(X)
M = ncol(Y)
N = nrow(Y)
nu= N-(M+F)+1 #nu must be positive
C_1 =  c(1  + x0 %*% xprimexinv %*% t(x0)) #for a prediction of the factor setting x0 (a vector of size F=ncol(X))
varY = A/(nu) 
postmean = x0 %*% hatB
nsim = 2000
ysim = rmvt(n=nsim,delta=postmux0,C_1*varY,df=nu) 

Теперь квантили ysim - это интервалы допуска бета-ожидания от прогнозирующего распределения, вы, конечно, можете напрямую использовать выборочное распределение, чтобы делать все, что вы хотите.

Пьер Лебрун
источник