<jupyter_start><jupyter_text>Modelamiento BayesianoEn el enfoque Bayesiano:- La inferencia se realiza calculado la distribución **a posteriori** de los parámetros- Se modela la incerteza de los datos y parámetros con una **distribución de probabilidad conjunta**- $\theta$ es un vector aleatorio cuya distribución de probabilidad se denomina **a priori**El teorema de Bayes y la ley de probabilidades totales nos permiten escribir:$$p(\theta| \cal{D}) = \frac{p(\cal{D}, \theta)}{p(\cal{D})}= \frac{p(\cal{D}|\theta) p(\theta)}{\int p(\cal{D}|\theta) p(\theta) d\theta} \propto p(\cal{D}|\theta) p(\theta),$$En el enfoque Bayesiano buscamos calcular la distribución a posteriori, la que se construye a partir de la función de verosimilitud, la distribución a priori y la evidencia, conocida también como la verosimilitud marginal de los datos.  El procedimiento de inferencia Bayesiana:1. Definir la función de verosimilitud de los datos2. Escoger una distribución a priori para los parámetros3. Construir la distribución conjunta de datos y parámetros4. Determinar la distribución a posteriori usando el Teorema de Bayes5. Calcular estimadores puntuales:**Moda a posteriori (MAP, maximum a posteriori)**$$\hat \theta = \text{arg} \max_\theta p(\theta|\cal{D}) =  \text{arg} \max_\theta p(\cal{D}| \theta) p(\theta),$$para esta maximización no se requiere la evidencia pues esta no depende de $\theta$. Además aplicando logaritmo podemos separar la verosimilitud de la distribución a priori:$$\hat \theta = \text{arg} \max_\theta (\log p(\cal{D}| \theta) + \log p(\theta)),$$**Esperanza a posteriori**$$\mathbb{E}(\theta \mid \cal{D}) = \int \theta p(\theta \mid \cal{D}) d\theta$$**Varianza a posteriori**$$\mathbb{V}(\theta \mid \cal{D}) = \int (\theta- \mathbb{E}(\theta \mid \cal{D}))^2 p(\theta \mid \cal{D}) d\theta$$  Algunas dificultades con MAP:- no es necesariamente representativa de la distribución a posteriori- no es invariante a reparametrizaciones de los parámetros<jupyter_code>library("IRdisplay")
display_png(file="graficoMAP0.png")
display_png(file="graficoMAP.png")<jupyter_output><empty_output><jupyter_text>Intervalos de credibilidad:- Intervalos centrales$$C_{\alpha}({\cal D}) = (l,u) : P(l \leq \theta \leq u \mid {\cal D}) = 1-\alpha$$- Regiones de mas altas densidades a posteriori (HPD):Sea $p*$ un umbral para la distribución a posteriori que cumple que:$$1-\alpha = \int_{\theta: p(\theta\mid {\cal D})>p*} p(\theta \mid {\cal D}) d\theta$$entonces se define:$$C_{\alpha}({\cal D}) = \{ \theta : p(\theta\mid {\cal D})\geq p*\}$$<jupyter_code>display_png(file="graficoCentral.png")<jupyter_output><empty_output><jupyter_text>**Distribución Predictiva a posteriori**$$p(x \mid \cal{D}) = \int p(x \mid \theta) p(\theta \mid \cal{D}) d\theta$$Nos permite predecir la probabilidad de ocurrencia de un nuevo dato $x$  Distribuciones a prioriDistribuciones naturales conjugadas: Una elección posible para las distribuciones a priori son aquellas que se conjugan con las funciones  de verosimilitud, de manera que la distribución a posteriori es del mismo tipo que la distribución a priori pero con nuevos parámetros. Veamos algunos ejemplos.**i) El modelo Beta-binomial**Consideremos el caso en que $$ X_i \sim i.i.d. \, Binomial(K,\theta) \qquad X_i \in \{0,...,K\}, i=1,\cdots,n$$$$ P(X_i = k\mid \theta) = {{K}\choose{k}} \theta^k (1-\theta)^{(K-k)} \qquad k = 1,...,K$$Y por otra parte$$ \theta \sim Beta(a,b) \qquad \theta \in [0,1], a,b>0$$$$f_{\theta}(\theta) = \frac{\Gamma(a+b)}{ \Gamma(a) \Gamma(b)} \theta^{a-1} (1-\theta)^{b-1} \qquad \theta \in [0,1]$$<jupyter_code>library(ggplot2)
options(repr.plot.width = 8, repr.plot.height = 6)
# Grafico densidad Beta con distintos parámetros
vec <- seq(0,1,by=0.01)
fvec <- data.frame(cbind(vec,dbeta(vec,0.1,0.1),dbeta(vec,1,1),dbeta(vec,2,3),dbeta(vec,8,4)))

h<- ggplot(fvec, aes(x=vec))
h <- h +coord_cartesian(ylim=c(0,3))
h <- h + geom_line(y=fvec$V5,color="green", linetype=1,size=0.8)
h <- h + geom_line(y=fvec$V4,color="red",size=0.8, linetype=2)
h <- h + geom_line(y=fvec$V3,color="blue",size=0.8, linetype=3)
h <- h + geom_line(y=fvec$V2,color="violet", size=0.8,linetype=4)
h<jupyter_output><empty_output><jupyter_text>**Algunas Propiedades de la Distribución Beta**i) $moda(\theta)= \underset{\theta}{\operatorname{argmax}} f_{\theta}(\theta)$$$\begin{align}\frac{\partial f_{\theta}(\theta)}{\partial \theta}=0 \qquad &\implies \qquad (a-1)\theta^{a-2}(1-\theta)^{b-1} -  (b-1)\theta^{a-1}(1-\theta)^{b-2} =0 \nonumber\\&\implies \qquad \theta^{a-2}(1-\theta)^{b-2}((a-1)(1-\theta) - (b-1)\theta) =0 \nonumber\\&\implies \qquad \theta(a-1 + b-1) =a-1 \nonumber\\&\implies \qquad moda(\theta) = \frac{a-1}{a+b-2} \nonumber\\\end{align}$$ii) $\mathbb{E}(\theta) = \int_0^1 \theta f_{\theta}(\theta) d\theta$$$\begin{align}\mathbb{E}(\theta) & = \frac{\Gamma(a+b)}{ \Gamma(a) \Gamma(b)} \int_0^1 \theta \theta^{a-1} (1-\theta)^{b-1}d\theta \nonumber\\& =  \frac{a}{a+b}\nonumber\\\end{align}$$iii) $\mathbb{V}(\theta) = \int_0^1 (\theta- \mathbb{E}(\theta)^2 f_{\theta}(\theta) d\theta$$$\begin{align}\mathbb{V}(\theta) & = \frac{\Gamma(a+b)}{ \Gamma(a) \Gamma(b)} \int_0^1 (\theta- \mathbb{E}(\theta)^2 \theta^{a-1} (1-\theta)^{b-1}d\theta \nonumber\\& =  \frac{ab}{(a+b)^2(a+b+1)}\nonumber\\\end{align}$$**Nota:**$$\Gamma(z) = \int_0^{\infty} t^{z-1} e^{-t} dt$$y cumple:$$\Gamma(z+1) = z \Gamma(z), \qquad \Gamma(n) = (n-1)! , n\in N$$ **Cálculo de la distribución a posteriori**$$\begin{align}f_{\theta}(\theta \mid \{X_i\}) &\propto \prod_{i=1}^n P(X_i=k_i\mid \theta)f_{\theta}(\theta) \nonumber \\&\propto \prod_{i=1}^n \theta^{k_i} (1-\theta)^{(K-k_i)} \theta^{a-1} (1-\theta)^{b-1} \nonumber \\&\nonumber\\&= \theta^{N_1} (1-\theta)^{N_0} \theta^{a-1} (1-\theta)^{b-1} \nonumber \\& \nonumber\\&= \theta^{N_1+a-1} (1-\theta)^{N_0+b-1}\end{align}$$donde $$N_1 = \sum_{i=1}^n k_i \qquad y \qquadN_0 = \sum_{i=1}^n (K-k_i)$$Es decir $$ \theta \mid \{X_i\}_{i=1,..,n} \sim Beta(N_1+a,N_0+b)$$**Moda a posteriori (MAP):**$$\theta_{MAP}= \underset{\theta}{\operatorname{argmax}} f_{\theta}(\theta \mid \{x_i\})= \frac{a+N_1-1}{a+b+N-2}$$**Media a posteriori:**$$\mathbb{E}(\theta \mid \{x_i\}_{i=1,...,n})  = \frac{a+N_1}{a+b+N}$$**Varianza a posteriori**$$\mathbb{V}(\theta \mid \{x_i\}_{i=1,...,n})  = \frac{(a+N_1)(b+N_2)}{(a+b+N)^2(a+b+N+1)}$$<jupyter_code>display_png(file="graficosBeta.png")<jupyter_output><empty_output><jupyter_text>**Ejercicio 1:** Calcule la distribución a posteriori para el caso del modelo Beta-Bernoulli.  **Ejercicio 2:** Calcule las siguientes probabilidades predictivas a posteriori:i) $P(x=1 \mid \cal{D})$ considerando que hay un sólo ensayo ( o sea K=1).ii)$P(x \mid \cal{D}, M)$ considerando que hay M ensayos. **i) El modelo Dirichlet-multinomial**Consideremos el caso en que $$ {\bf{X}}=(x_1,....x_K) \sim Multinomial(n,{\bf{\theta}}) \qquad x_i \in \{0,...,n\}, n=\sum_{k=1}^K x_k , \qquad {\bf{\theta}}=(\theta_1,...,\theta_K)  $$$$ P(X = (x_1,...,x_k) \mid {\bf{\theta} }) = {{n}\choose{x_1,...,x_k}}\prod_{k=1}^K\theta_k^{x_k} = \frac{n!}{x_1!...x_k!}\prod_{k=1}^K\theta_k^{x_k}$$De manera que la función de verosimilitud para n observaciones queda:$$ P(\cal{D} \mid \bf{\theta}) \propto \prod_{k=1}^K\theta_k^{N_k}$$donde$$N_k = \sum_{i=1}^n (X_i)_k$$Por otra parte, sea$$ \theta = (\theta_1,...\theta_K) \sim Dirichlet(\alpha_1,...,\alpha_K) \qquad \theta_i \in [0,1], \sum_{k=1}^K \theta_i= 1, \alpha_i>0$$$$f_{\theta}(\theta) = \frac{1}{B(\alpha)} \prod_{k=1}^K \theta_k^{\alpha_k-1} \qquad B(\alpha) = \frac{\prod_{k=1}^K \Gamma(\alpha_k)}{\Gamma(\alpha_1+...+\alpha_K)}$$En este caso, resulta que la distribución a posteriori es también Dirichlet:$$\begin{align}f_{\theta}(\theta \mid \cal{D}) & \propto P(\cal{D} \mid \theta)f_{\theta}(\theta) \nonumber \\& \propto \prod_{k=1}^K \theta_k^{N_k}\theta_k^{\alpha_k-1} = \prod_{k=1}^K \theta_k^{N_k+\alpha_k-1} \nonumber \end{align}$$Es decir$$(\theta \mid {\cal{D}}) \sim Dirichlet(\alpha_1+N_1,...,\alpha_K+N_k)$$  Ejemplo: diferencia de proporcionesSupongamos que estamos interesados en comprar en una plataforma online, y hay dos vendedores que entregan el producto al mismo precio. El vendedor 1 tiene 90 evaluaciones positivas y 10 negativas. El vendendor 2 tiene 2 evaluaciones positivas y 0 negativas. ¿Cómo decidir cual comprar?En el enfoque bayesiano, el análisis es como sigue:Sean $\theta_1$ y $\theta_2$ las confiabilidades reales de los dos vendedores. Como no tenemos mas información, vamos a considerar a priori uniformes: $\theta_i \sim Beta(1,1)$ y entonces podemos usar el modelo beta-binomial. De manera que$$\theta_i \mid {\cal D} \sim Beta(1+y_i, 1 + N_i - y_i), \qquad i=1,2$$En este caso $$N_1 = 100, \,y_1 = 90,\qquad N_2 = 2,\, y_2 =2$$Entonces podemos calcular, para $\delta = \theta_1 - \theta_2$:$$P(\delta > 0 \mid {\cal D}) = \int_0^1 \int_0^1 \mathbb{1}_{\{\theta_1 - \theta_2>0\} }Beta(1+y_1, 1 + N_1 - y_1)$$$$\qquad \qquad \qquad \qquad Beta(1+y_2, 1 + N_2 - y_2)d\theta_1 d\theta_2$$ Utilizando integración numérica se puede obtener que $P(\delta > 0 \mid {\cal D})= 0.71$ y un intervalo central de la figura.<jupyter_code>display_png(file="proporciones.png")<jupyter_output><empty_output><jupyter_text>Enfoque Bayesiano para la Selección de modelosSi disponemos de un conjunto de modelos entre los cuales elegir, por ejemplo una familia de distribuciones de probabilidad paramétricas de complejidad diferente, estamos en presencia del problema de selección de modelos.<jupyter_code>display_png(file="overfitting.png")<jupyter_output><empty_output><jupyter_text>En el enfoque frecuentista, una manera de abordar esto es usar validación cruzada (cross-validation), para definir los errores para cada modelo y luego seleccionar el de mínimo error. En el enfoque Bayesiano podemos hacer uso de la distribución predictiva de los datos o evidencia, por cada modelo $m$ y escoger el de mayor probabilidad dados los datos $\cal{D}$$$p(m \mid \cal{D}) = \frac{p(\cal{D} \mid m) p(m) }{\sum_{m \in \cal{M}}p(\cal{D} \mid m) p(m)}$$y entonces calculamos el modelo MAP:$$\hat{m} = \underset{m \in \cal{M}}{\operatorname{argmax}}p(m \mid \cal{D})$$Si consideramos una distribución  a priori no informativa sobre los modelos (distribución uniforme), es decir $p(m) \propto 1$, tenemos que:$$p(m \mid \cal{D}) \propto p(\cal{D} \mid m) = \int p(\cal{D} \mid \theta) p(\theta \mid m) d\theta$$llamada verosimilitud marginal o evidencia del modelo $m$. **La navaja de Occam:** (Bayesian Occam's razor): el hecho de que la verosimilitud marginal considere la integración respecto de los parámetros impide que se maximice simplemente por aumentar el número de parámetros, con lo cual se previene el sobreajuste en la elección del modelo.<jupyter_code>display_png(file="seleccion1.png")
display_png(file="seleccion2.png")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.198176627682107, 'fraction_numerical': 0.020720139976056727, 'mean_word_length': 4.957213384530993, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '352357', 'license': 'Unlicense', 'n_tokens_mistral': 4337, 'n_tokens_neox': 4183, 'n_words': 1069, 'path': '10_ModelamientoBayesiano/ModelamientoBayesiano.ipynb', 'repo_name': 'magister-informatica-uach/INFO337', 'size': '10859'}