Bewerten Sie die posteriore prädiktive Verteilung in der Bayes'schen linearen Regression

Ich bin verwirrt darüber, wie die posteriore prädiktive Verteilung für die Bayes'sche lineare Regression nach dem hier auf Seite 3 beschriebenen und unten kopierten Grundfall bewertet werden soll.

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

Der Grundfall ist dieses lineare Regressionsmodell:

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

Wenn wir entweder einen einheitlichen Prior auf mit einem Scale-Inv vor oder den Normal-Inverse-Gamma-Prior (siehe hier ) verwenden, ist die posteriore Vorhersageverteilung analytisch und ist Student t. χ 2 σ $\beta$$\chi^2$$\sigma^2$

Was ist mit diesem Modell?

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

Wenn , aber bekannt ist, ist die posteriore Vorhersageverteilung multivariate Gaußsche. Normalerweise kennen Sie , müssen es aber schätzen. Vielleicht sagst du seine Diagonale und machst die Diagonale in irgendeiner Weise zu einer Funktion der Kovariaten. Dies wird im Kapitel über die lineare Regression der Bayes'schen Datenanalyse von Gelman erörtert .Σ $y \sim N(X\beta, \Sigma)$$\Sigma$$\Sigma$

Gibt es in diesem Fall eine analytische Form für die posteriore Vorhersageverteilung? Kann ich meine Schätzung einfach in einen multivariaten Schüler t einfügen? Wenn Sie mehr als eine Varianz schätzen, ist die Verteilung immer noch multivariate Schüler t?

Ich frage , weil ich einige sagen haben bereits auf der Hand. Ich möchte wissen, ob es wahrscheinlicher ist, dass es beispielsweise durch lineare Regression A, lineare Regression B vorhergesagt wurde $\tilde y$

Wenn Sie eine einheitliche vor auf nehmen , dann die hintere für ist mit Um die prädiktive Verteilung zu finden, benötigen wir weitere Informationen. Wenn und bedingt unabhängig von bei , dann ist Aber normalerweise sind und für diese nicht bedingt unabhängig, sondern haben wir normalerweise β β | y ~ N ( β , V β ) . Β = [ X ' Σ - 1 X ] X '$\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).$$ Wenn dies der Fall ist, dann Dies setzt voraus, dass und bekannt sind. Wie Sie hervorheben, sind sie normalerweise unbekannt und müssen geschätzt werden. Für die gängigen Modelle mit dieser Struktur, z. B. Zeitreihen und räumliche Modelle, gibt es im Allgemeinen keine geschlossene Form für die Vorhersageverteilung. $$\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}$

Unter nicht informativen oder multivariaten Normal-Wishart-Prioritäten haben Sie die analytische Form als multivariate Schülerverteilung für eine klassische mutlivariate multiple Regression. Ich denke, die Entwicklungen in diesem Dokument hängen mit Ihrer Frage zusammen (vielleicht gefällt Ihnen Anhang A :-)). Normalerweise habe ich das Ergebnis mit einer posterioren Vorhersageverteilung verglichen, die mit WinBUGS und der analytischen Form erhalten wurde: Sie sind genau gleichwertig. Das Problem wird nur dann schwierig, wenn Sie in Modellen mit gemischten Effekten zusätzliche zufällige Effekte haben, insbesondere bei unausgeglichenem Design.

Im Allgemeinen sind bei klassischen Regressionen y und ỹ bedingt unabhängig (Residuen sind iid)! Wenn dies nicht der Fall ist, ist die hier vorgeschlagene Lösung natürlich nicht korrekt.

In R (hier Lösung für einheitliche Prioritäten) wird unter der Annahme, dass Sie ein lm-Modell (mit dem Namen "Modell") einer der Antworten in Ihrem Modell erstellt und als "Modell" bezeichnet haben, wie Sie die multivariate Vorhersageverteilung erhalten

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) 

Jetzt sind Quantile von ysim Beta-Erwartungstoleranzintervalle von der Vorhersageverteilung. Sie können die Stichprobenverteilung natürlich direkt verwenden, um zu tun, was Sie wollen.