Regresión de Bayes: ¿cómo se hace en comparación con la regresión estándar?

Question

Tengo algunas preguntas sobre la regresión Bayesiana:

1. Dada una regresión estándar como $y = \beta_0 + \beta_1 x + \varepsilon$. Si quiero convertir esto en una regresión Bayesiana, ¿necesito distribuciones previas tanto para $\beta_0$ como para $\beta_1$ (o no funciona de esta manera)?

2. En la regresión estándar uno intentaría minimizar los residuos para obtener valores únicos para $\beta_0$ y $\beta_1. ¿Cómo se hace esto en la regresión de Bayes?

---

Realmente estoy luchando mucho aquí:

$$ \text{posterior} = \text{prior} \times \text{likelihood} $$

La probabilidad proviene del conjunto de datos actual (entonces es mi parámetro de regresión pero no como un único valor sino como una distribución de probabilidad, ¿verdad?). La priori proviene de una investigación previa (digamos). Así que tengo esta ecuación:

$$ y = \beta_1 x + \varepsilon $$

con $\beta_1$ siendo mi probabilidad o posterior (¿o esto es completamente erróneo)?

Simplemente no puedo entender cómo se transforma la regresión estándar en una de Bayes.

Answer 1

El modelo de regresión lineal simple

$$ y_i = \alpha + \beta x_i + \varepsilon $$

puede ser escrito en términos del modelo probabilístico que lo respalda

$$ \mu_i = \alpha + \beta x_i \\ y_i \sim \mathcal{N}(\mu_i, \sigma) $$

es decir la variable dependiente $Y$ sigue una distribución normal parametrizada por la media $\mu_i$, que es una función lineal de $X$ parametrizada por $\alpha,\beta$, y por la desviación estándar $\sigma$. Si estimas un modelo así utilizando mínimos cuadrados ordinarios, no necesitas preocuparte por la formulación probabilística, porque estás buscando los valores óptimos de los parámetros $\alpha,\beta$ minimizando los errores cuadrados entre los valores ajustados y los valores predichos. Por otro lado, podrías estimar dicho modelo utilizando estimación de máxima verosimilitud, donde buscarías los valores óptimos de los parámetros maximizando la función de verosimilitud

$$ \DeclareMathOperator*{\argmax}{arg\,max} \argmax_{\alpha,\,\beta,\,\sigma} \prod_{i=1}^n \mathcal{N}(y_i; \alpha + \beta x_i, \sigma) $$

donde $\mathcal{N}$ es una función de densidad de la distribución normal evaluada en los puntos de $y_i$, parametrizada por las medias $\alpha + \beta x_i$ y la desviación estándar $\sigma$.

En el enfoque bayesiano en lugar de maximizar la función de verosimilitud sola, asumiríamos distribuciones previas para los parámetros y usaríamos el teorema de Bayes

$$ \text{posterior} \propto \text{likelihood} \times \text{prior} $$

La función de verosimilitud es la misma que arriba, pero lo que cambia es que asumes algunas distribuciones previas para los parámetros estimados $\alpha,\beta,\sigma$ y las incluyes en la ecuación

$$ \underbrace{f(\alpha,\beta,\sigma\mid Y,X)}_{\text{posterior}} \propto \underbrace{\prod_{i=1}^n \mathcal{N}(y_i\mid \alpha + \beta x_i, \sigma)}_{\text{likelihood}} \; \underbrace{f_{\alpha}(\alpha) \, f_{\beta}(\beta) \, f_{\sigma}(\sigma)}_{\text{priors}} $$

"¿Qué distribuciones?" es una pregunta diferente, ya que hay un número ilimitado de opciones. Para los parámetros $\alpha,\beta$ podrías, por ejemplo, asumir distribuciones normales parametrizadas por algunos hiperparámetros, o una distribución $t$ si deseas asumir colas más pesadas, o una distribución uniforme si no deseas hacer demasiadas suposiciones, pero quieres asumir que los parámetros pueden ser "cualquier cosa en un rango dado" a priori, etc. Para $\sigma$ necesitas asumir alguna distribución prior que esté limitada a ser mayor que cero, ya que la desviación estándar debe ser positiva. Esto puede llevar a la formulación del modelo como se ilustra a continuación por John K. Kruschke.

(fuente: http://www.indiana.edu/~kruschke/BMLR/)

Mientras que en la máxima verosimilitud estabas buscando un único valor óptimo para cada uno de los parámetros, en el enfoque bayesiano aplicando el teorema de Bayes obtienes la distribución posterior de los parámetros. La estimación final dependerá de la información que proviene de tus datos y de tus priors, pero cuanta más información contenga tus datos, menos influyentes son los priors.

Observa que al usar priors uniformes, toman la forma $f(\theta) \propto 1$ después de eliminar las constantes de normalización. Esto hace que el teorema de Bayes sea proporcional solo a la función de verosimilitud, por lo que la distribución posterior alcanzará su máximo exactamente en el mismo punto que la estimación de máxima verosimilitud. Por lo tanto, la estimación bajo priors uniformes será la misma que al utilizar mínimos cuadrados ordinarios ya que minimizar los errores cuadrados corresponde a maximizar la verosimilitud normal.

Para estimar un modelo en el enfoque bayesiano en algunos casos puedes usar priors conjugados, por lo que la distribución posterior está directamente disponible (ver ejemplo aquí). Sin embargo, en la gran mayoría de los casos, la distribución posterior no estará directamente disponible y tendrás que usar métodos de Monte Carlo de cadena de Markov para estimar el modelo (consulta este ejemplo de cómo usar el algoritmo Metropolis-Hastings para estimar parámetros de regresión lineal). Finalmente, si solo estás interesado en estimaciones puntuales de los parámetros, podrías usar estimación del máximo a posteriori, es decir

$$ \argmax_{\alpha,\,\beta,\,\sigma} f(\alpha,\beta,\sigma\mid Y,X) $$

Para una descripción más detallada de la regresión logística, puedes consultar el hilo Modelo de logit bayesiano: ¿explicación intuitiva?.

Para aprender más, puedes consultar los siguientes libros:

Kruschke, J. (2014). Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan. Academic Press.

Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2004). Bayesian data analysis. Chapman & Hall/CRC.

Answer 2

Dado un conjunto de datos $D = (x_1,y_1), \ldots, (x_N,y_N)$ donde $x \in \mathbb{R}^d, y \in \mathbb{R}$, un modelo de Regresión Lineal Bayesiano plantea el problema de la siguiente manera:

Prior: $$w \sim \mathcal{N}(0, \sigma_w^2 I_d)$$

$w$ es un vector $(w_1, \ldots, w_d)^T$, por lo que la distribución anterior es una Gaussiana multivariante; y $I_d$ es la matriz identidad $d\times d$.

Verosimilitud: $$Y_i \sim \mathcal{N}(w^T x_i, \sigma^2)$$

Asumimos que $Y_i \perp Y_j | w, i \neq j$

Por ahora usaremos la precisión en lugar de la varianza, $a = 1/\sigma^2$, y $b = 1/\sigma_w^2$. También asumiremos que $a,b$ son conocidos.

El prior puede ser expresado como $$p(w) \propto \exp \Big\{ -\frac{b}{2} w^t w \Big\}$$

Y la verosimilitud $$p(D|w) \propto \exp \Big\{ -\frac{a}{2} (y-Aw)^T (y-Aw) \Big\}$$

donde $y = (y_1,\ldots,y_N)^T$ y $A$ es una matriz $n\times d$ donde la fila i-ésima es $x_i^T$.

Entonces el posterior es $$p(w|D) \propto p(D|w) p(w)$$

Después de muchos cálculos descubrimos que

$$p(w|D) \sim \mathcal{N}(w | \mu, \Lambda^{-1})$$

donde ($\Lambda$ es la matriz de precisión)

$$\Lambda = a A^T A + b I_d $$ $$\mu = a \Lambda^{-1} A^T y$$

Observa que $\mu$ es igual a $w_{MAP}$ de la regresión lineal regular, esto es porque para la Distribución Gaussiana, la media es igual a la moda.

También, podemos hacer algunas operaciones algebraicas en $\mu$ y obtener la siguiente igualdad ($\Lambda = aA^TA+bI_d$):

$$\mu = (A^T A + \frac{b}{a} I_d)^{-1} A^T y$$

y comparar con $w_{MLE}$:

$$w_{MLE} = (A^T A)^{-1} A^T y$$

La expresión adicional en $\mu$ corresponde al prior. Esto es similar a la expresión para la regresión Ridge, para el caso especial cuando $\lambda = \frac{b}{a}$. La regresión Ridge es más general porque la técnica puede elegir priors impropios (en la perspectiva Bayesiana).

Para la distribución predictiva posterior:

$$p(y|x,D) = \int p(y|x,D,w) p(w|x,D) dw = \int p(y|x,w) p(w|D) dw$$

es posible calcular que

$$y|x,D \sim \mathcal{N}(\mu^Tx, \frac{1}{a} + x^T \Lambda^{-1}x)$$

Referencia: Lunn et al. El libro de BUGS

Para usar una herramienta MCMC como JAGS/Stan revisa Doing Bayesian Data Analysis de Kruschke