Oceń tylny rozkład predykcyjny w regresji liniowej Bayesa

Jestem zdezorientowany, jak ocenić tylny rozkład predykcyjny dla regresji liniowej Bayesa, pomijając podstawowy przypadek opisany tutaj na stronie 3, i skopiowałem poniżej.

$$p(\tilde y \mid y) = \int p(\tilde y \mid \beta, \sigma^2) p(\beta, \sigma^2 \mid y)$$

Podstawowym przypadkiem jest ten model regresji liniowej:

$$y = X \beta + \epsilon, \hspace{10mm} y \sim N(X \beta, \sigma^2)$$

Jeśli użyjemy albo munduru przed , ze skalą Inv przed , LUB normalnego-odwrotnego-gamma przed (patrz tutaj ), tylna dystrybucja predykcyjna jest analityczna i jest studentem. $\beta$$\chi^2$$\sigma^2$

Co z tym modelem?

$$y = X \beta + \epsilon, \hspace{10mm} y \sim N(X \beta, \Sigma)$$

Kiedy , ale jest znane, tylny rozkład predykcyjny jest wielowymiarowy Gaussa. Zwykle nie znasz , ale musisz to oszacować. Może powiesz, że jest przekątna, i uczynisz ją w jakiś sposób funkcją zmiennych towarzyszących. Jest to omówione w rozdziale dotyczącym regresji liniowej w analizie danych bayesowskich Gelmana .$y \sim N(X\beta, \Sigma)$$\Sigma$$\Sigma$

Czy w tym przypadku istnieje forma analityczna dla tylnego rozkładu predykcyjnego? Czy mogę po prostu podłączyć swoją ocenę tego do wielowymiarowego ucznia t? Jeśli oszacujesz więcej niż jedną wariancję, czy rozkład nadal jest wielowymiarowy dla studenta t?

Pytam, bo powiedz, że mam już trochę . Chcę wiedzieć, czy jest bardziej prawdopodobne, że można to przewidzieć np. Przez regresję liniową A, regresję liniową B. $\tilde y$

Jeśli założysz mundur przed , to tylna dla to with Aby znaleźć rozkład predykcyjny, potrzebujemy więcej informacji. Jeśli i jest warunkowo niezależny od dającego , to Ale zazwyczaj dla tego typu modeli i nie są warunkowo niezależne, zamiast tego zazwyczaj mamy $\beta$$\beta$

$$\beta|y \sim N(\hat{\beta},V_\beta).$$ $$\hat{\beta} = [X'\Sigma^{-1}X]X'y \quad \mbox{and} \quad V_\beta =[X'\Sigma^{-1}X]^{-1}.$$ $\tilde{y} \sim N(\tilde{X}\beta,\tilde{\Sigma})$$y$$\beta$ $$\tilde{y}|y \sim N(\tilde{X}\hat{\beta}, \tilde{\Sigma} + V_\beta).$$ $y$$\tilde{y}$ $$\left( \begin{array}{c} y \\ \tilde{y} \end{array} \right) \sim N\left( \left[\begin{array}{c} X\beta \\ \tilde{X}\beta \end{array}\right], \left[ \begin{array}{cc} \Sigma & \Sigma_{12} \\ \Sigma_{21} & \tilde{\Sigma} \end{array} \right]\right).$$ W takim przypadku Zakłada się, że i są znane. Jak zauważyłeś, zazwyczaj są one nieznane i należy je oszacować. W przypadku popularnych modeli o tej strukturze, np. Szeregów czasowych i modeli przestrzennych, na ogół nie będzie zamkniętej formy rozkładu predykcyjnego. $$\tilde{y}|y \sim N(\tilde{X}\hat{\beta} + \Sigma_{21}\Sigma^{-1}(y-X\hat{\beta}), \tilde{\Sigma} - \Sigma_{21}\Sigma^{-1}\Sigma_{12}).$$ $\Sigma, \Sigma_{12},$$\tilde{\Sigma}$

W przypadku nieinformacyjnych lub wielowymiarowych priorów normalnego życzenia, masz formę analityczną jako wielowymiarowy rozkład Studenta, dla klasycznej wielowariantowej, regresji wielokrotnej. Myślę, że zmiany w tym dokumencie są związane z twoim pytaniem (może ci się spodobać załącznik A :-)). Zazwyczaj porównywałem wynik z późniejszym rozkładem predykcyjnym uzyskanym za pomocą WinBUGS i formą analityczną: są one dokładnie równoważne. Problem staje się trudny tylko wtedy, gdy masz dodatkowe efekty losowe w modelach z efektami mieszanymi, szczególnie w przypadku niezrównoważonego projektu.

Ogólnie, przy klasycznych regresjach, yi ỹ są warunkowo niezależne (reszty są iid)! Oczywiście, jeśli tak nie jest, proponowane tutaj rozwiązanie jest nieprawidłowe.

W R (tutaj rozwiązanie dla jednolitych priorytetów), zakładając, że wykonałeś model lm (o nazwie „model”) jednej z odpowiedzi w swoim modelu i nazwałeś go „modelem”, oto jak uzyskać wielowymiarowy rozkład predykcyjny

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) 

Teraz kwantyle ysim są przedziałami tolerancji oczekiwań beta od rozkładu predykcyjnego, możesz oczywiście bezpośrednio użyć próbkowanego rozkładu, aby zrobić co chcesz.