Regresión lineal multivariante bayesiana

Enfoque bayesiano para la regresión lineal multivariante

En estadística , la regresión lineal multivariante bayesiana es un enfoque bayesiano de la regresión lineal multivariante , es decir, una regresión lineal en la que el resultado previsto es un vector de variables aleatorias correlacionadas en lugar de una única variable aleatoria escalar. Se puede encontrar un tratamiento más general de este enfoque en el artículo Estimador MMSE .

Detalles

Consideremos un problema de regresión donde la variable dependiente que se va a predecir no es un único escalar de valor real sino un vector de longitud m de números reales correlacionados. Como en la configuración de regresión estándar, hay n observaciones, donde cada observación i consta de k −1 variables explicativas , agrupadas en un vector de longitud k (donde se ha añadido una variable ficticia con un valor de 1 para permitir un coeficiente de intersección). Esto puede verse como un conjunto de m problemas de regresión relacionados para cada observación i : donde el conjunto de errores están todos correlacionados. De manera equivalente, puede verse como un único problema de regresión donde el resultado es un vector de fila y los vectores de coeficientes de regresión se apilan uno al lado del otro, de la siguiente manera: ${\displaystyle \mathbf {x} _{i}}$ $${\displaystyle {\begin{aligned}y_{i,1}&=\mathbf {x} _{i}^{\mathsf {T}}{\boldsymbol {\beta }}_{1}+\epsilon _{i,1}\\&\;\;\vdots \\y_{i,m}&=\mathbf {x} _{i}^{\mathsf {T}}{\boldsymbol {\beta }}_{m}+\epsilon _{i,m}\end{aligned}}}$$ ${\displaystyle \{\epsilon _{i,1},\ldots ,\epsilon _{i,m}\}}$ ${\displaystyle \mathbf {y} _{i}^{\mathsf {T}}}$ $${\displaystyle \mathbf {y} _{i}^{\mathsf {T}}=\mathbf {x} _{i}^{\mathsf {T}}\mathbf {B} +{\boldsymbol {\epsilon }}_{i}^{\mathsf {T}}.}$$

La matriz de coeficientes B es una matriz donde los vectores de coeficientes para cada problema de regresión se apilan horizontalmente: ${\displaystyle k\times m}$ ${\displaystyle {\boldsymbol {\beta }}_{1},\ldots ,{\boldsymbol {\beta }}_{m}}$ $${\displaystyle \mathbf {B} ={\begin{bmatrix}{\begin{pmatrix}\\{\boldsymbol {\beta }}_{1}\\\\\end{pmatrix}}\cdots {\begin{pmatrix}\\{\boldsymbol {\beta }}_{m}\\\\\end{pmatrix}}\end{bmatrix}}={\begin{bmatrix}{\begin{pmatrix}\beta _{1,1}\\\vdots \\\beta _{k,1}\end{pmatrix}}\cdots {\begin{pmatrix}\beta _{1,m}\\\vdots \\\beta _{k,m}\end{pmatrix}}\end{bmatrix}}.}$$

El vector de ruido para cada observación i es conjuntamente normal, de modo que los resultados para una observación dada están correlacionados: ${\displaystyle {\boldsymbol {\epsilon }}_{i}}$ $${\displaystyle {\boldsymbol {\epsilon }}_{i}\sim N(0,{\boldsymbol {\Sigma }}_{\epsilon }).}$$

Podemos escribir todo el problema de regresión en forma matricial como: donde Y y E son matrices. La matriz de diseño X es una matriz con las observaciones apiladas verticalmente, como en la configuración de regresión lineal estándar : $${\displaystyle \mathbf {Y} =\mathbf {X} \mathbf {B} +\mathbf {E} ,}$$ ${\displaystyle n\times m}$ ${\displaystyle n\times k}$ $${\displaystyle \mathbf {X} ={\begin{bmatrix}\mathbf {x} _{1}^{\mathsf {T}}\\\mathbf {x} _{2}^{\mathsf {T}}\\\vdots \\\mathbf {x} _{n}^{\mathsf {T}}\end{bmatrix}}={\begin{bmatrix}x_{1,1}&\cdots &x_{1,k}\\x_{2,1}&\cdots &x_{2,k}\\\vdots &\ddots &\vdots \\x_{n,1}&\cdots &x_{n,k}\end{bmatrix}}.}$$

La solución clásica, frecuentista, de mínimos cuadrados lineales consiste simplemente en estimar la matriz de coeficientes de regresión utilizando la pseudoinversa de Moore-Penrose : ${\displaystyle {\hat {\mathbf {B} }}}$ $${\displaystyle {\hat {\mathbf {B} }}=(\mathbf {X} ^{\mathsf {T}}\mathbf {X} )^{-1}\mathbf {X} ^{\mathsf {T}}\mathbf {Y} .}$$

Para obtener la solución bayesiana, necesitamos especificar la probabilidad condicional y luego encontrar la probabilidad conjugada previa adecuada. Al igual que con el caso univariante de la regresión bayesiana lineal , descubriremos que podemos especificar una probabilidad conjugada previa condicional natural (que depende de la escala).

Escribamos nuestra probabilidad condicional como ^[1] escribiendo el error en términos de y obtenemos $${\displaystyle \rho (\mathbf {E} |{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp \left(-{\tfrac {1}{2}}\operatorname {tr} \left(\mathbf {E} ^{\mathsf {T}}\mathbf {E} {\boldsymbol {\Sigma }}_{\epsilon }^{-1}\right)\right),}$$ ${\displaystyle \mathbf {E} }$ ${\displaystyle \mathbf {Y} ,\mathbf {X} ,}$ ${\displaystyle \mathbf {B} }$ $${\displaystyle \rho (\mathbf {Y} |\mathbf {X} ,\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {Y} -\mathbf {X} \mathbf {B} )^{\mathsf {T}}(\mathbf {Y} -\mathbf {X} \mathbf {B} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1})),}$$

Buscamos una conjugación natural previa, es decir, una densidad conjunta que tenga la misma forma funcional que la probabilidad. Como la probabilidad es cuadrática en , reescribimos la probabilidad de modo que sea normal en (la desviación de la estimación de la muestra clásica). ${\displaystyle \rho (\mathbf {B} ,\Sigma _{\epsilon })}$ ${\displaystyle \mathbf {B} }$ ${\displaystyle (\mathbf {B} -{\hat {\mathbf {B} }})}$

Utilizando la misma técnica que en la regresión lineal bayesiana , descomponemos el término exponencial utilizando una forma matricial de la técnica de suma de cuadrados. Sin embargo, aquí también necesitaremos utilizar el cálculo diferencial de matrices ( producto de Kronecker y transformaciones de vectorización ).

Primero, apliquemos la suma de cuadrados para obtener una nueva expresión para la probabilidad: $${\displaystyle \rho (\mathbf {Y} |\mathbf {X} ,\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-(n-k)/2}\exp(-\operatorname {tr} ({\tfrac {1}{2}}\mathbf {S} ^{\mathsf {T}}\mathbf {S} {\boldsymbol {\Sigma }}_{\epsilon }^{-1}))|{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {B} -{\hat {\mathbf {B} }})^{\mathsf {T}}\mathbf {X} ^{\mathsf {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})),}$$ $${\displaystyle \mathbf {S} =\mathbf {Y} -\mathbf {X} {\hat {\mathbf {B} }}}$$

Nos gustaría desarrollar una forma condicional para las probabilidades anteriores: donde es una distribución Wishart inversa y es alguna forma de distribución normal en la matriz . Esto se logra utilizando la transformación de vectorización , que convierte la probabilidad de una función de las matrices a una función de los vectores . $${\displaystyle \rho (\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })=\rho ({\boldsymbol {\Sigma }}_{\epsilon })\rho (\mathbf {B} |{\boldsymbol {\Sigma }}_{\epsilon }),}$$ ${\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon })}$ ${\displaystyle \rho (\mathbf {B} |{\boldsymbol {\Sigma }}_{\epsilon })}$ ${\displaystyle \mathbf {B} }$ ${\displaystyle \mathbf {B} ,{\hat {\mathbf {B} }}}$ ${\displaystyle {\boldsymbol {\beta }}=\operatorname {vec} (\mathbf {B} ),{\hat {\boldsymbol {\beta }}}=\operatorname {vec} ({\hat {\mathbf {B} }})}$

Escribir $${\displaystyle \operatorname {tr} ((\mathbf {B} -{\hat {\mathbf {B} }})^{\mathsf {T}}\mathbf {X} ^{\mathsf {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})=\operatorname {vec} (\mathbf {B} -{\hat {\mathbf {B} }})^{\mathsf {T}}\operatorname {vec} (\mathbf {X} ^{\mathsf {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})}$$

Sea donde denota el producto de Kronecker de las matrices A y B , una generalización del producto externo que multiplica una matriz por una matriz para generar una matriz, que consiste en cada combinación de productos de elementos de las dos matrices. $${\displaystyle \operatorname {vec} (\mathbf {X} ^{\mathsf {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})=({\boldsymbol {\Sigma }}_{\epsilon }^{-1}\otimes \mathbf {X} ^{\mathsf {T}}\mathbf {X} )\operatorname {vec} (\mathbf {B} -{\hat {\mathbf {B} }}),}$$ ${\displaystyle \mathbf {A} \otimes \mathbf {B} }$ ${\displaystyle m\times n}$ ${\displaystyle p\times q}$ ${\displaystyle mp\times nq}$

Entonces esto conducirá a una probabilidad que es normal en . $${\displaystyle {\begin{aligned}&\operatorname {vec} (\mathbf {B} -{\hat {\mathbf {B} }})^{\mathsf {T}}({\boldsymbol {\Sigma }}_{\epsilon }^{-1}\otimes \mathbf {X} ^{\mathsf {T}}\mathbf {X} )\operatorname {vec} (\mathbf {B} -{\hat {\mathbf {B} }})\\&=({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\mathsf {T}}({\boldsymbol {\Sigma }}_{\epsilon }^{-1}\otimes \mathbf {X} ^{\mathsf {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})\end{aligned}}}$$ ${\displaystyle ({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})}$

Con la probabilidad en una forma más manejable, ahora podemos encontrar un conjugado previo natural (condicional).

Distribución previa conjugada

El conjugado natural anterior que utiliza la variable vectorizada tiene la forma: ^[1] donde y ${\displaystyle {\boldsymbol {\beta }}}$ $${\displaystyle \rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon })=\rho ({\boldsymbol {\Sigma }}_{\epsilon })\rho ({\boldsymbol {\beta }}|{\boldsymbol {\Sigma }}_{\epsilon }),}$$ $${\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon })\sim {\mathcal {W}}^{-1}(\mathbf {V} _{0},{\boldsymbol {\nu }}_{0})}$$ $${\displaystyle \rho ({\boldsymbol {\beta }}|{\boldsymbol {\Sigma }}_{\epsilon })\sim N({\boldsymbol {\beta }}_{0},{\boldsymbol {\Sigma }}_{\epsilon }\otimes {\boldsymbol {\Lambda }}_{0}^{-1}).}$$

Distribución posterior

Utilizando la distribución previa y de probabilidad anteriores, la distribución posterior se puede expresar como: ^[1] donde . Los términos que involucran se pueden agrupar (con ) utilizando: con $${\displaystyle {\begin{aligned}\rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\propto {}&|{\boldsymbol {\Sigma }}_{\epsilon }|^{-({\boldsymbol {\nu }}_{0}+m+1)/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} (\mathbf {V} _{0}{\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}\\&\times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {B} -\mathbf {B} _{0})^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} -\mathbf {B} _{0}){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}\\&\times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {Y} -\mathbf {XB} )^{\mathsf {T}}(\mathbf {Y} -\mathbf {XB} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))},\end{aligned}}}$$ ${\displaystyle \operatorname {vec} (\mathbf {B} _{0})={\boldsymbol {\beta }}_{0}}$ ${\displaystyle \mathbf {B} }$ ${\displaystyle {\boldsymbol {\Lambda }}_{0}=\mathbf {U} ^{\mathsf {T}}\mathbf {U} }$ $${\displaystyle {\begin{aligned}&\left(\mathbf {B} -\mathbf {B} _{0}\right)^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}\left(\mathbf {B} -\mathbf {B} _{0}\right)+\left(\mathbf {Y} -\mathbf {XB} \right)^{\mathsf {T}}\left(\mathbf {Y} -\mathbf {XB} \right)\\={}&\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {U} \mathbf {B} _{0}\end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} \right)^{\mathsf {T}}\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {U} \mathbf {B} _{0}\end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} \right)\\={}&\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {U} \mathbf {B} _{0}\end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} _{n}\right)^{\mathsf {T}}\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {U} \mathbf {B} _{0}\end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} _{n}\right)+\left(\mathbf {B} -\mathbf {B} _{n}\right)^{\mathsf {T}}\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}\right)\left(\mathbf {B} -\mathbf {B} _{n}\right)\\={}&\left(\mathbf {Y} -\mathbf {X} \mathbf {B} _{n}\right)^{\mathsf {T}}\left(\mathbf {Y} -\mathbf {X} \mathbf {B} _{n}\right)+\left(\mathbf {B} _{0}-\mathbf {B} _{n}\right)^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}\left(\mathbf {B} _{0}-\mathbf {B} _{n}\right)+\left(\mathbf {B} -\mathbf {B} _{n}\right)^{\mathsf {T}}\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}\right)\left(\mathbf {B} -\mathbf {B} _{n}\right),\end{aligned}}}$$ $${\displaystyle \mathbf {B} _{n}=\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}\right)^{-1}\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} {\hat {\mathbf {B} }}+{\boldsymbol {\Lambda }}_{0}\mathbf {B} _{0}\right)=\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}\right)^{-1}\left(\mathbf {X} ^{\mathsf {T}}\mathbf {Y} +{\boldsymbol {\Lambda }}_{0}\mathbf {B} _{0}\right).}$$

Esto ahora nos permite escribir la posterior en una forma más útil: $${\displaystyle {\begin{aligned}\rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\propto {}&|{\boldsymbol {\Sigma }}_{\epsilon }|^{-({\boldsymbol {\nu }}_{0}+m+n+1)/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {V} _{0}+(\mathbf {Y} -\mathbf {XB_{n}} )^{\mathsf {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B} _{n}-\mathbf {B} _{0})^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} _{n}-\mathbf {B} _{0})){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}\\&\times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {B} -\mathbf {B} _{n})^{\mathsf {T}}(\mathbf {X} ^{T}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})(\mathbf {B} -\mathbf {B} _{n}){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}.\end{aligned}}}$$

Esto toma la forma de una distribución Wishart inversa multiplicada por una distribución normal matricial : y $${\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\sim {\mathcal {W}}^{-1}(\mathbf {V} _{n},{\boldsymbol {\nu }}_{n})}$$ $${\displaystyle \rho (\mathbf {B} |\mathbf {Y} ,\mathbf {X} ,{\boldsymbol {\Sigma }}_{\epsilon })\sim {\mathcal {MN}}_{k,m}(\mathbf {B} _{n},{\boldsymbol {\Lambda }}_{n}^{-1},{\boldsymbol {\Sigma }}_{\epsilon }).}$$

Los parámetros de esta posterior están dados por: $${\displaystyle \mathbf {V} _{n}=\mathbf {V} _{0}+(\mathbf {Y} -\mathbf {XB_{n}} )^{\mathsf {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B} _{n}-\mathbf {B} _{0})^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} _{n}-\mathbf {B} _{0})}$$ $${\displaystyle {\boldsymbol {\nu }}_{n}={\boldsymbol {\nu }}_{0}+n}$$ $${\displaystyle \mathbf {B} _{n}=(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\mathsf {T}}\mathbf {Y} +{\boldsymbol {\Lambda }}_{0}\mathbf {B} _{0})}$$ $${\displaystyle {\boldsymbol {\Lambda }}_{n}=\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}}$$

Véase también

• Regresión lineal bayesiana
• Distribución normal matricial

Referencias

1. Peter E. Rossi, Greg M. Allenby, Rob McCulloch. Estadística bayesiana y marketing . John Wiley & Sons, 2012, pág. 32.

• Box, GEP ; Tiao, GC (1973). "8". Inferencia bayesiana en análisis estadístico . Wiley. ISBN 0-471-57428-7.
• Geisser, S. (1965). "Estimación bayesiana en análisis multivariante". Anales de estadística matemática . 36 (1): 150–159. JSTOR  2238083.
• Tiao, GC; Zellner, A. (1964). "Sobre la estimación bayesiana de la regresión multivariante". Revista de la Royal Statistical Society. Serie B (Metodológica) . 26 (2): 277–285. JSTOR  2984424.