Évaluer la distribution prédictive postérieure dans la régression linéaire bayésienne

Je ne sais pas comment évaluer la distribution prédictive postérieure de la régression linéaire bayésienne, après le cas de base décrit ici à la page 3, et copié ci-dessous.

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

Le cas de base est ce modèle de régression linéaire:

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

Si nous utilisons soit un a priori uniforme sur , avec une échelle Inv- priori sur , OU le a priori gamma normal-inverse (voir ici ), la distribution prédictive postérieure est analytique et est t de Student. $\beta$$\chi^2$$\sigma^2$

Et pour ce modèle?

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

Lorsque , mais est connu, la distribution prédictive postérieure est gaussienne multivariée. Habituellement, vous ne connaissez pas , mais devez l'estimer. Vous dites peut-être sa diagonale et faites en sorte que la diagonale soit fonction des covariables d'une manière ou d'une autre. Ceci est discuté dans le chapitre sur la régression linéaire de l'analyse des données bayésiennes de Gelman .$y \sim N(X\beta, \Sigma)$$\Sigma$$\Sigma$

Existe-t-il une forme analytique pour la distribution prédictive postérieure dans ce cas? Puis-je simplement brancher mon estimation de celui-ci dans un étudiant multivarié t? Si vous estimez plus d'une variance, la distribution est-elle toujours multivariée t étudiant?

Je demande parce que disons que j'ai déjà sous la main. Je veux savoir si elle est plus susceptible d'avoir été prédite par exemple par régression linéaire A, régression linéaire B $\tilde y$

Si vous supposez un précédent uniforme sur , alors le postérieur de est avec Pour trouver la distribution prédictive, nous avons besoin de plus d'informations. Si et est conditionnellement indépendant de donné , alors Mais typiquement pour ces types de modèles, et ne sont pas conditionnellement indépendants, à la place, nous avons généralement $\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).$$ Si tel est le cas, alors Cela suppose que et sont tous connus. Comme vous le faites remarquer, ils sont généralement inconnus et doivent être estimés. Pour les modèles communs qui ont cette structure, par exemple les séries chronologiques et les modèles spatiaux, il n'y aura généralement pas de forme fermée pour la distribution prédictive. $$\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}$

Sous les a priori Normal-Wishart non informatifs ou multivariés, vous avez la forme analytique en tant que distribution de Student multivariée, pour une régression multiple mutlivariée classique. Je suppose que les développements dans ce document sont liés à votre question (vous aimerez peut-être l'annexe A :-)). J'ai généralement comparé le résultat avec une distribution prédictive postérieure obtenue en utilisant WinBUGS et la forme analytique: ils sont exactement équivalents. Le problème ne devient difficile que lorsque vous avez des effets aléatoires supplémentaires dans les modèles à effets mixtes, en particulier dans la conception asymétrique.

En général, avec les régressions classiques, y et ỹ sont conditionnellement indépendants (les résidus sont iid)! Bien sûr, si ce n'est pas le cas, la solution proposée ici n'est pas correcte.

Dans R, (ici, solution pour des a priori uniformes), en supposant que vous avez créé un modèle lm (nommé "modèle") de l'une des réponses dans votre modèle, et que vous l'avez appelé "modèle", voici comment obtenir la distribution prédictive multivariée

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) 

Maintenant, les quantiles de ysim sont des intervalles de tolérance d'espérance bêta de la distribution prédictive, vous pouvez bien sûr utiliser directement la distribution échantillonnée pour faire ce que vous voulez.