$A\rightarrow H$: Una hipotesis (por ejemplo: Los parametros de un modelo son $\hat{\theta}$)
$B\rightarrow E$: Eventos/Evidencia/Observaciones
$$ P(H\vert E) = \frac{P(E\vert H)P(H)}{P(E)} $$$P(E)$ Se calcula con el teorema de probabilidades totales:
$$ P(E) = \sum_i P(E\vert H_i)P(H_i) $$$P(E)$ Es la parte mas dificil de calcular, pero corresponde a una normalizacion, asi que si se puede evaluar $P(E\vert H)$ y $P(H) \implies$ Se puede integrar y normalizar.
$P(H\vert E) = \dfrac{P(E\vert H)P(H)}{P(E)}$
Nota: Gaussiana multivariada (Utiles para definir los priors).
En general: $$ G\left(\vec{X}\right) = \dfrac{1}{\sqrt{(2\pi)^n\text{det}(\Sigma)}}\times e^{\left(-\frac{1}{2}\left(\vec{x}-\vec{\mu}\right)^T\Sigma^{-1}\left(\vec{x}-\vec{\mu}\right)\right)} $$
Note que para el caso $N=1\rightarrow\Sigma=\sigma^2$ $$ \implies G(X) = \dfrac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} $$
$(x_i, y_i)_{i=1\ldots N}$ son las observaciones
$y_i = \beta_0 + \beta_1x_i + \epsilon_i \leftarrow$ modelo lineal
Se define
$$ \vec{y} = \left(\begin{matrix}y_1 \\ \vdots \\ y_N \end{matrix}\right) \qquad \vec{\epsilon} = \left(\begin{matrix}\epsilon_1 \\ \vdots \\ \epsilon_N \end{matrix}\right) \qquad \vec{\beta} = \left(\begin{matrix}\beta_1 \\ \beta_1 \end{matrix}\right) $$Con $\vec{\beta}$ los parametros del modelo.
Asumiendo $\epsilon_i = \epsilon = const$
$$\begin{matrix} P(\vec{x},\vec{y}\vert\vec{\beta}) = \prod_{i=1}^N\frac{1}{\sqrt{2\pi\epsilon^2}}e^{-\frac{\left(y_i - \beta_0 - \beta_1x_i\right)^2}{2\epsilon^2}} \\ \\ = \frac{1}{(2\pi\epsilon^2)^{N/2}}e^{-\frac{1}{2\epsilon^2}\sum_{i=1}^N\left(y_i - \beta_0 - \beta_1x_i\right)^2} \\ \\ P\left(\vec{X},\vec{Y}\vert\vec{\beta}\right) = \frac{1}{(2\pi\epsilon^2)^{N/2}}e^{-\frac{1}{2\epsilon^2}\left(\vec{y}-\bar{X}\vec{\beta}\right)^T\left(\vec{y}-\bar{X}\vec{\beta}\right)}\end{matrix} $$Si el prior es gaussiano y la verosimilitud tambien $ \implies P(E\vert H)P(H) $ tambien es gaussiana, luego se simplifica el problema utilizando propiedades analiticas conocidas.
In [1]:
import scipy as sp
import matplotlib.pyplot as plt
%matplotlib notebook
In [2]:
x, y = sp.mgrid[-3:3:99j, -3:4:99j]
In [3]:
def Gauss2D(x, y, mu, sigma_pars):
mu_x, mu_y = my
sigma_x, sigma_y, rho = sigma_pars
A = 1. / (2. * sp.pi * 1. * 1. * sp.sqrt(1. - .5**2))
B = sp.exp(-((x - mu_x) ** 2 / sigma_x**2 + (y - mu_y) ** 2) / sigma_y**2 -
2 * rho * (x - mu_x) * (y - mu_y) / sigma_x / sigma_y) / (2 * (1 - rho**2))
return A * B