Evaluar la distribución predictiva posterior en regresión lineal bayesiana

10

Estoy confundido sobre cómo evaluar la distribución predictiva posterior para la regresión lineal bayesiana, más allá del caso básico descrito aquí en la página 3, y copiado a continuación.

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

El caso básico es este modelo de regresión lineal:

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

Si usamos un previo uniforme en , con una escala-Inv anterior en , O el anterior normal-inverso-gamma (ver aquí ) la distribución predictiva posterior es analítica y es t de Student. χ 2 σ 2βχ2σ2

¿Qué hay de este modelo?

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

Cuando , pero se conoce, la distribución predictiva posterior es gaussiana multivariante. Por lo general, no sabes , pero tienes que estimarlo. Tal vez dices que es diagonal y haces que la diagonal sea una función de las covariables de alguna manera. Esto se discute en el capítulo de regresión lineal del Análisis de datos bayesianos de Gelman .Σ ΣyN(Xβ,Σ)ΣΣ

¿Existe una forma analítica para la distribución predictiva posterior en este caso? ¿Puedo conectar mi estimación en una t de estudiante multivariante? Si estima más de una varianza, ¿la distribución sigue siendo multivariante t de estudiante?

Me pregunto porque dicen que tengo un poco de ya en la mano. Quiero saber si es más probable que haya sido predicho por ejemplo, por regresión lineal A, regresión lineal B y~

bill_e
fuente
1
Si tiene muestras posteriores de la distribución posterior, puede evaluar la distribución predictiva con una aproximación de Monte Carlo
niandra82
Ah, gracias, siempre podría hacer eso. ¿No hay una fórmula analítica en este caso?
bill_e
Los enlaces están rotos, por cierto. Sería genial si incorporaras las referencias de otra manera.
Maxim.K

Respuestas:

6

Si asume un uniforme anterior en , entonces el posterior para es with Para encontrar la distribución predictiva, necesitamos más información. Si y es condicionalmente independiente de dado , entonces Pero, por lo general, para este tipo de modelos, y no son condicionalmente independientes, en su lugar, normalmente tenemos β β | y ~ N ( β , V β ) . Β = [ X ' Σ - 1 X ] X ' yββ

β|yN(β^,Vβ).
β^=[XΣ1X]XyandVβ=[XΣ1X]1.
y~N(X~β,Σ~)yβ
y~|yN(X~β^,Σ~+Vβ).
yy~˜y| y~N(~Xβ+Σ21Σ-1(Y-Xβ),~Σ-Σ21Σ-1Σ12). Σ,Σ12,˜Σ
(yy~)N([XβX~β],[ΣΣ12Σ21Σ~]).
Si este es el caso, entonces Esto supone que y son todos conocidos. Como señala, generalmente son desconocidos y deben estimarse. Para los modelos comunes que tienen esta estructura, por ejemplo, series de tiempo y modelos espaciales, generalmente no habrá una forma cerrada para la distribución predictiva.
y~|yN(X~β^+Σ21Σ1(yXβ^),Σ~Σ21Σ1Σ12).
Σ,Σ12,Σ~
jaradniemi
fuente
2

En las versiones previas no informativas o multivariadas de Normal-Wishart, tiene la forma analítica como una distribución de estudiante multivariada, para una regresión múltiple, clásica, multivariada. Supongo que los desarrollos en este documento están relacionados con su pregunta (puede que le guste el Apéndice A :-)). Típicamente comparé el resultado con una distribución predictiva posterior obtenida usando WinBUGS y la forma analítica: son exactamente equivalentes. El problema solo se vuelve difícil cuando tiene efectos aleatorios adicionales en modelos de efectos mixtos, especialmente en diseños desequilibrados.

¡En general, con regresiones clásicas, y y ỹ son condicionalmente independientes (los residuos son iid)! Por supuesto, si no es el caso, entonces la solución propuesta aquí no es correcta.

En R, (aquí, la solución para los antecedentes uniformes), suponiendo que usted hizo un modelo lm (llamado "modelo") de una de las respuestas en su modelo, y lo llamó "modelo", aquí es cómo obtener la distribución predictiva multivariada

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) 

Ahora, los cuantiles de ysim son intervalos de tolerancia de expectativa beta de la distribución predictiva, por supuesto, puede usar directamente la distribución muestreada para hacer lo que quiera.

Pierre Lebrun
fuente