Regresión lineal bayesiana

Compartir Imprimir Citar

La regresión lineal bayesiana es un tipo de modelado condicional en el que la media de una variable se describe mediante una combinación lineal de otras variables, con el objetivo de obtener la probabilidad posterior de los coeficientes de regresión (así como otros parámetros que describen la distribución de la regresora).) y, en última instancia, permitir la predicción fuera de la muestra del regresor (a menudo etiquetado como y) condicional a los valores observados de los regresores (generalmente X). La versión más simple y más utilizada de este modelo es el modelo lineal normal, en el que ydadoXse distribuye Gaussiana. En este modelo, y bajo una elección particular de probabilidades previas para los parámetros, los llamados anteriores conjugados, el posterior se puede encontrar analíticamente. Con priores elegidos más arbitrariamente, los posteriores generalmente tienen que ser aproximados.

Configuración del modelo

Considere un problema de regresión lineal estándar, en el i=1,ldots,nque especificamos la media de la distribución condicional de y_{yo}un k  veces 1vector predictor dado mathbf {x} _{i}:

{displaystyle y_{i}=mathbf {x} _{i}^{mathsf {T}}{boldsymbol {beta }}+varepsilon _{i},}

donde { símbolo de negrita { beta}}es un k  veces 1vector, y varepsilon _{i}son variables aleatorias independientes e idénticamente distribuidas normalmente:

{displaystyle varepsilon _{i}sim N(0,sigma ^{2}).}

Esto corresponde a la siguiente función de verosimilitud:

{displaystyle rho (mathbf {y} mid mathbf {X},{boldsymbol {beta }},sigma ^{2})propto (sigma ^{2})^{-n/ 2}exp left(-{frac {1}{2sigma ^{2}}}(mathbf {y} -mathbf {X} {boldsymbol {beta }})^{mathsf { T}}(mathbf {y} -mathbf {X} {boldsymbol {beta }})right).}

La solución de mínimos cuadrados ordinarios se usa para estimar el vector de coeficientes usando el pseudoinverso de Moore-Penrose:

{displaystyle {hat {boldsymbol {beta }}}=(mathbf {X} ^{mathsf {T}}mathbf {X})^{-1}mathbf {X} ^{mathsf {T}}matemáticas {y} }

donde mathbf{X}está la n veces kmatriz de diseño, cada fila de la cual es un vector predictor {displaystyle mathbf {x}_{i}^{mathsf {T}}}; y mathbf {y}es el nortevector- columna {displaystyle [y_{1};cdots ;y_{n}]^{mathsf {T}}}.

Este es un enfoque frecuentista y asume que hay suficientes medidas para decir algo significativo acerca de { símbolo de negrita { beta}}. En el enfoque bayesiano, los datos se complementan con información adicional en forma de una distribución de probabilidad previa. La creencia previa sobre los parámetros se combina con la función de verosimilitud de los datos según el teorema de Bayes para producir la creencia posterior sobre los parámetros { símbolo de negrita { beta}}y sigma. El previo puede tomar diferentes formas funcionales dependiendo del dominio y la información que está disponible a priori.

Dado que los datos comprenden tanto mathbf {y}como mathbf{X}, el enfoque solo en la distribución de mathbf {y}condicional en la mathbf{X}justificación de necesidades. De hecho, un análisis bayesiano "completo" requeriría una probabilidad conjunta {displaystyle rho (mathbf {y},mathbf {X} mid {boldsymbol {beta }},sigma ^{2},gamma)}junto con un anterior { estilo de visualización  rho ( beta,  sigma ^ {2},  gamma)}, donde gamasimboliza los parámetros de la distribución para mathbf{X}. Solo bajo el supuesto de exogeneidad (débil) se puede factorizar la probabilidad conjunta en {displaystyle rho (mathbf {y} mid {boldsymbol {mathbf {X} }},beta,sigma ^{2})rho (mathbf {X} mid gamma)}. La última parte generalmente se ignora bajo el supuesto de conjuntos de parámetros disjuntos. Más aún, bajo los supuestos clásicos mathbf{X}se considera elegido (por ejemplo, en un experimento diseñado) y por lo tanto tiene una probabilidad conocida sin parámetros.

Con antecedentes conjugados

Distribución previa conjugada

Para una distribución previa arbitraria, puede que no haya una solución analítica para la distribución posterior. En esta sección, consideraremos un denominado conjugado anterior para el cual la distribución posterior se puede derivar analíticamente.

Un prior rho(boldsymbolbeta,sigma^{2})es conjugado a esta función de verosimilitud si tiene la misma forma funcional con respecto a { símbolo de negrita { beta}}y sigma. Dado que el logaritmo de la verosimilitud es cuadrático en { símbolo de negrita { beta}}, el logaritmo de la verosimilitud se vuelve a escribir de tal manera que la probabilidad se vuelve normal en (boldsymbolbeta-sombrero{boldsymbolbeta}). Escribe

{displaystyle {begin{alineado}(mathbf {y} -mathbf {X} {boldsymbol {beta }})^{mathsf {T}}(mathbf {y} -mathbf {X} {boldsymbol {beta }})&=[(mathbf {y} -mathbf {X} {hat {boldsymbol {beta }}})+(mathbf {X} {hat {boldsymbol {beta }}}-mathbf {X} {boldsymbol {beta }})]^{mathsf {T}}[(mathbf {y} -mathbf {X} {hat {boldsymbol { beta }}})+(mathbf {X} {hat {boldsymbol {beta }}}-mathbf {X} {boldsymbol {beta }})]\&=(mathbf {y} -mathbf {X} {hat {boldsymbol {beta }}})^{mathsf {T}}(mathbf {y} -mathbf {X} {hat { boldsymbol {beta }}})+({boldsymbol {beta }}-{hat {boldsymbol {beta }}})^{mathsf {T}}(mathbf {X} ^{ mathsf {T}}mathbf {X})({boldsymbol {beta }}-{hat {boldsymbol {beta }}})+underbrace {2(mathbf {X} {hat { boldsymbol {beta }}}-mathbf {X} {boldsymbol {beta }})^{mathsf {T}}(mathbf {y} -mathbf {X} {hat {boldsymbol { beta }}})} _ {= 0}\&=(mathbf {y} -mathbf {X} {hat {boldsymbol {beta }}})^{mathsf {T}}(mathbf {y} -mathbf {X} {hat {boldsymbol {beta }}})+({boldsymbol {beta }}-{hat {boldsymbol {beta }}})^{ mathsf {T}}(mathbf {X} ^{mathsf {T}}mathbf {X})({boldsymbol {beta }}-{hat {boldsymbol {beta }}}),.end{alineado}}}=(mathbf {y} -mathbf {X} {hat {boldsymbol {beta }}})^{mathsf {T}}(mathbf {y} -mathbf {X} {hat { boldsymbol {beta }}})+({boldsymbol {beta }}-{hat {boldsymbol {beta }}})^{mathsf {T}}(mathbf {X} ^{ mathsf {T}}mathbf {X})({boldsymbol {beta }}-{hat {boldsymbol {beta }}}),.end{alineado}}}=(mathbf {y} -mathbf {X} {hat {boldsymbol {beta }}})^{mathsf {T}}(mathbf {y} -mathbf {X} {hat { boldsymbol {beta }}})+({boldsymbol {beta }}-{hat {boldsymbol {beta }}})^{mathsf {T}}(mathbf {X} ^{ mathsf {T}}mathbf {X})({boldsymbol {beta }}-{hat {boldsymbol {beta }}}),.end{alineado}}}

La probabilidad ahora se reescribe como

{displaystyle rho (mathbf {y} |mathbf {X},{boldsymbol {beta }},sigma ^{2})propto (sigma ^{2})^{-{frac {v}{2}}}exp left(-{frac {vs^{2}}{2{sigma }^{2}}}right)(sigma ^{2})^{- {frac {nv}{2}}}exp left(-{frac {1}{2{sigma }^{2}}}({boldsymbol {beta }}-{hat { boldsymbol {beta }}})^{mathsf {T}}(mathbf {X} ^{mathsf {T}}mathbf {X})({boldsymbol {beta }}-{hat { boldsymbol {beta }}})derecha),}

dónde

{displaystyle vs^{2}=(mathbf {y} -mathbf {X} {hat {boldsymbol {beta }}})^{mathsf {T}}(mathbf {y} - mathbf {X} {hat {boldsymbol {beta }}})quad {text{ y }}quad v=nk,}

donde kes el número de coeficientes de regresión.

Esto sugiere una forma para el anterior:

{displaystyle rho ({boldsymbol {beta }},sigma ^{2})=rho (sigma ^{2})rho ({boldsymbol {beta }}mid sigma ^{ 2}),}

donde { estilo de visualización  rho ( sigma ^ {2})}es una distribución gamma inversa

{displaystyle rho (sigma ^{2})propto (sigma ^{2})^{-{frac {v_{0}}{2}}-1}exp left(-{ fracción {v_{0}s_{0}^{2}}{2sigma ^{2}}}right).}

En la notación introducida en el artículo de distribución de gamma inversa, esta es la densidad de una {displaystyle {text{Inv-Gamma}}(a_{0},b_{0})}distribución con {displaystyle a_{0}={tfrac {v_{0}}{2}}}y {displaystyle b_{0}={tfrac{1}{2}}v_{0}s_{0}^{2}}con v_{0}y {displaystyle s_{0}^{2}}como los valores previos de vy ^{2}, respectivamente. De manera equivalente, también se puede describir como una distribución de chi-cuadrado inversa escalada,{displaystyle {text{Escala-inv-}}chi ^{2}(v_{0},s_{0}^{2}).}

Además, la densidad previa condicional rho(boldsymbolbeta|sigma^{2})es una distribución normal,

{displaystyle rho ({boldsymbol {beta }}mid sigma ^{2})propto (sigma ^{2})^{-k/2}exp left(-{frac { 1}{2sigma ^{2}}}({boldsymbol {beta }}-{boldsymbol {mu }}_{0})^{mathsf {T}}mathbf {Lambda } _ {0}({boldsymbol {beta }}-{boldsymbol {mu }}_{0})right).}

En la notación de la distribución normal, la distribución previa condicional es{displaystyle {mathcal {N}}left({boldsymbol {mu }}_{0},sigma ^{2}{boldsymbol {Lambda }}_{0}^{-1} Correcto).}

Distribución posterior

Con el ahora anterior especificado, la distribución posterior se puede expresar como

{displaystyle {begin{alineado}rho ({boldsymbol {beta }},sigma ^{2}mid mathbf {y},mathbf {X})&propto rho (mathbf { y} mid mathbf {X},{boldsymbol {beta }},sigma ^{2})rho ({boldsymbol {beta }}mid sigma ^{2})rho ( sigma ^{2})\&propto (sigma ^{2})^{-n/2}exp left(-{frac {1}{2{sigma }^{2}}} (mathbf {y} -mathbf {X} {boldsymbol {beta }})^{mathsf {T}}(mathbf {y} -mathbf {X} {boldsymbol {beta }}) right)(sigma ^{2})^{-k/2}exp left(-{frac {1}{2sigma ^{2}}}({boldsymbol {beta }}- {boldsymbol {mu}}_{0})^{mathsf {T}}{boldsymbol {Lambda}}_{0}({boldsymbol {beta}}-{boldsymbol {mu} }_{0})right)(sigma ^{2})^{-(a_{0}+1)}exp left(-{frac {b_{0}}{sigma ^{2 }}}right)end{alineado}}}

Con alguna reorganización, el posterior se puede reescribir para que la media posterior negritamu_ndel vector de parámetros { símbolo de negrita { beta}}se pueda expresar en términos del estimador de mínimos cuadrados hat{boldsymbolbeta}y la media anterior { símbolo de negrita { mu}}_{0}, con la fuerza de la anterior indicada por la matriz de precisión previa.boldsymbolLambda_0

{displaystyle {boldsymbol {mu }}_{n}=(mathbf {X} ^{mathsf {T}}mathbf {X} +{boldsymbol {Lambda }}_{0})^ {-1}(mathbf {X} ^{mathsf {T}}mathbf {X} {hat {boldsymbol {beta }}}+{boldsymbol {Lambda }}_{0}{ símbolo de negrita {mu }}_{0}).}

Para justificar que negritamu_nde hecho es la media posterior, los términos cuadráticos en la exponencial se pueden reorganizar como una forma cuadrática en {displaystyle {boldsymbol {beta }}-{boldsymbol {mu }}_{n}}.

{displaystyle (mathbf {y} -mathbf {X} {boldsymbol {beta }})^{mathsf {T}}(mathbf {y} -mathbf {X} {boldsymbol {beta }})+({boldsymbol {beta }}-{boldsymbol {mu }}_{0})^{mathsf {T}}{boldsymbol {Lambda }}_{0}({ boldsymbol {beta }}-{boldsymbol {mu }}_{0})=({boldsymbol {beta }}-{boldsymbol {mu }}_{n})^{mathsf {T }}(mathbf {X} ^{mathsf {T}}mathbf {X} +{boldsymbol {Lambda }}_{0})({boldsymbol {beta }}-{boldsymbol { mu }}_{n})+mathbf {y} ^{mathsf {T}}mathbf {y} -{boldsymbol {mu }}_{n}^{mathsf {T}}( mathbf {X} ^{mathsf {T}}mathbf {X} +{boldsymbol {Lambda }}_{0}){boldsymbol {mu }}_{n}+{boldsymbol {mu }}_{0}^{mathsf {T}}{boldsymbol {Lambda }}_{0}{boldsymbol {mu }}_{0}.}

Ahora, el posterior se puede expresar como una distribución normal multiplicada por una distribución gamma inversa:

{displaystyle rho ({boldsymbol {beta }},sigma ^{2}mid mathbf {y},mathbf {X})propto (sigma ^{2})^{-k/ 2}exp left(-{frac {1}{2{sigma }^{2}}}({boldsymbol {beta }}-{boldsymbol {mu }}_{n})^ {mathsf {T}}(mathbf {X} ^{mathsf {T}}mathbf {X} +mathbf {Lambda }_{0})({boldsymbol {beta }}-{ símbolo de negrita {mu }}_{n})right)(sigma ^{2})^{-{frac {n+2a_{0}}{2}}-1}exp left(-{ frac {2b_{0}+mathbf {y} ^{mathsf {T}}mathbf {y} -{boldsymbol {mu }}_{n}^{mathsf {T}}(mathbf {X} ^{mathsf {T}}mathbf {X} +{boldsymbol {Lambda }}_{0}){boldsymbol {mu }}_{n}+{boldsymbol {mu } }_{0}^{mathsf {T}}{boldsymbol {Lambda }}_{0}{boldsymbol {mu }}_{0}}{2sigma ^{2}}}right).}

Por lo tanto, la distribución posterior se puede parametrizar de la siguiente manera.

{displaystyle rho ({boldsymbol {beta }},sigma ^{2}mid mathbf {y},mathbf {X})propto rho ({boldsymbol {beta }}mid sigma ^{2},mathbf {y},mathbf {X})rho (sigma ^{2}mid mathbf {y},mathbf {X}),}

donde los dos factores corresponden a las densidades mathcal{N}left(boldsymbolmu_n, sigma^2boldsymbolLambda_n^{-1} right),y text{Inv-Gamma}left(a_n,b_n right)distribuciones, con los parámetros de estos dados por

{displaystyle {boldsymbol {Lambda }}_{n}=(mathbf {X} ^{mathsf {T}}mathbf {X} +mathbf {Lambda }_{0}),quad {boldsymbol {mu }}_{n}=({boldsymbol {Lambda }}_{n})^{-1}(mathbf {X} ^{mathsf {T}}mathbf {X } {sombrero {boldsymbol {beta }}}+{boldsymbol {Lambda }}_{0}{boldsymbol {mu }}_{0}),}
{displaystyle a_{n}=a_{0}+{frac {n}{2}},qquad b_{n}=b_{0}+{frac {1}{2}}(mathbf { y} ^{mathsf {T}}mathbf {y} +{boldsymbol {mu }}_{0}^{mathsf {T}}{boldsymbol {Lambda }}_{0}{ boldsymbol {mu }}_{0}-{boldsymbol {mu }}_{n}^{mathsf {T}}{boldsymbol {Lambda }}_{n}{boldsymbol {mu } }_{norte}).}

Esto se puede interpretar como un aprendizaje bayesiano donde los parámetros se actualizan de acuerdo con las siguientes ecuaciones.

{displaystyle {boldsymbol {mu }}_{n}=left(mathbf {X} ^{mathsf {T}}mathbf {X} +{boldsymbol {Lambda }}_{0} right)^{-1}left({boldsymbol {Lambda }}_{0}{boldsymbol {mu }}_{0}+mathbf {X} ^{mathsf {T}} mathbf {X} { sombrero { símbolo de negrita { beta}}}  derecho),}
{displaystyle {boldsymbol {Lambda }}_{n}=(mathbf {X} ^{mathsf {T}}mathbf {X} +{boldsymbol {Lambda }}_{0}), }
{displaystyle a_{n}=a_{0}+{frac {n}{2}},}
{displaystyle b_{n}=b_{0}+{frac {1}{2}}(mathbf {y} ^{mathsf {T}}mathbf {y} +{boldsymbol {mu } }_{0}^{mathsf {T}}{boldsymbol {Lambda }}_{0}{boldsymbol {mu }}_{0}-{boldsymbol {mu }}_{n} ^{mathsf {T}}{boldsymbol {Lambda }}_{n}{boldsymbol {mu }}_{n}).}

Evidencia modelo

La evidencia del modelo {displaystyle p(mathbf {y} mid m)}es la probabilidad de los datos dado el modelo metro. También se conoce como probabilidad marginal y como densidad predictiva previa. Aquí, el modelo está definido por la función de verosimilitud {displaystyle p(mathbf {y} mid mathbf {X},{boldsymbol {beta }},sigma)}y la distribución previa de los parámetros, es decirp(símbolo en negritabeta,sigma). La evidencia del modelo captura en un solo número qué tan bien dicho modelo explica las observaciones. La evidencia del modelo de regresión lineal bayesiana presentada en esta sección se puede utilizar para comparar modelos lineales competidores mediante la comparación de modelos bayesianos. Estos modelos pueden diferir en el número y valores de las variables predictoras, así como en sus antecedentes sobre los parámetros del modelo. La evidencia del modelo ya tiene en cuenta la complejidad del modelo, porque margina los parámetros integrando {displaystyle p(mathbf {y},{boldsymbol {beta }},sigma mid mathbf {X})}todos los valores posibles de { símbolo de negrita { beta}}y sigma.

{displaystyle p(mathbf {y} |m)=int p(mathbf {y} mid mathbf {X},{boldsymbol {beta }},sigma),p({boldsymbol {beta}},sigma),d{boldsymbol {beta}},dsigma}

Esta integral se puede calcular analíticamente y la solución se da en la siguiente ecuación.

{displaystyle p(mathbf {y} mid m)={frac {1}{(2pi)^{n/2}}}{sqrt {frac {det({boldsymbol { Lambda }}_{0})}{det({boldsymbol {Lambda }}_{n})}}}cdot {frac {b_{0}^{a_{0}}}{b_{ n}^{a_{n}}}}cdot {frac {Gamma (a_{n})}{Gamma (a_{0})}}}

Aquí Gamadenota la función gamma. Debido a que hemos elegido un conjugado anterior, la probabilidad marginal también se puede calcular fácilmente evaluando la siguiente igualdad para valores arbitrarios de { símbolo de negrita { beta}}y sigma.

{displaystyle p(mathbf {y} mid m)={frac {p({boldsymbol {beta }},sigma |m),p(mathbf {y} mid mathbf {X },{boldsymbol {beta }},sigma,m)}{p({boldsymbol {beta }},sigma mid mathbf {y},mathbf {X},m)}}}

Tenga en cuenta que esta ecuación no es más que una reorganización del teorema de Bayes. Insertando las fórmulas para la anterior, la verosimilitud y la posterior y simplificando la expresión resultante se obtiene la expresión analítica dada anteriormente.

Otros casos

En general, puede ser imposible o poco práctico derivar analíticamente la distribución posterior. Sin embargo, es posible aproximar el posterior mediante un método de inferencia bayesiano aproximado, como el muestreo de Monte Carlo o el Bayes variacional.

El caso especial boldsymbolmu_0=0, mathbf{Lambda}_0 = cmathbf{I}se llama regresión de cresta.

Se puede realizar un análisis similar para el caso general de la regresión multivariante y parte de esto proporciona la estimación bayesiana de matrices de covarianza: consulte Regresión lineal multivariante bayesiana.