# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats

# # Ley de los grandes números para variables no iid
#
# Previamente vimos que el promedio de un conjunto de $N$ variables independientes e idénticamente distribuidas (iid) converge a su valor esperado cuando $N$ es grande
#
# $$
# \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^N X_i = \mu
# $$
#
# También vimos que la cadena de markov, un proceso estocástico donde no se cumple el supuesto iid, puede tener en ciertos casos una distribución estacionaria
#
# **Recordatorio:** La distribución estacionaria $\pi$ de una cadena de Markov con matriz de transición $P$ es tal que $\pi P = \pi$
#
# ## Teorema de ergodicidad
#
# Una cadena de Markov irreducible y aperiodica tiene una distribución estacionaria $\pi$ única, independiente de valor del estado inicial y que cumple
#
# $$
# \lim_{n\to \infty} s_j(n) = \pi_j
# $$
#
# donde los componentes de $\pi$ representan la fracción de tiempo que la cadena estará en cada uno de los estados luego de observarla por un largo tiempo
#
# El límite de observar la cadena por un tiempo largo es análogo al análisis de estadísticos estáticos sobre muestras grandes. Esto es el equivalente a la ley de los grandes números para el caso de la cadena de Markov
#
#
# ## Nota histórica
#
# <NAME> mostró la primera versión de la Ley de los grandes números en su Ars Conjectandi en 1713. Esta primera versión parte del supuesto de que las VAs son iid. Bernoulli era un firme creyente del destino, se oponía al libre albedrío y abogaba por el determinismo en los fenómenos aleatorios
#
# En 1913 el matemático ruso <NAME> celebró el bicentenario de la famosa prueba de Bernoulli organizando un simposio donde presentó su nueva versión de la Ley de los grandes números que aplica sobre la clase de procesos estocásticos que hoy llamamos procesos de Markov, de esta forma extendiendo el resultado de Bernoulli a un caso que no es iid
#
# ## Más sobre la pugna de Markov y Nekrasov
#
# En aquellos tiempos Markov estaba en pugna con otro matemático ruso: <NAME>. Nekrasov había publicado previamente que "la independencia era una condición necesaria para la ley de los grandes números" y que la ley de los grandes números, que se observa en estadísticas sociales, asegura entonces que las personas actuan voluntariamente y con libre albedrío. Markov reaccionó a esta afirmación desarrollando un contra-ejemplo que terminó siendo lo que hoy conocemos como los procesos de Markov 
#
#

# # Principales usos de las cadenas de Markov
#
# Las cadenas de Markov tienen dos usos principales
#
# En primer lugar las cadenas de Markov se ocupan como **modelo o aproximación de fenómenos que evolucionan en el tiempo**. Esto es lo que vimos la lección anterior.
#
# En estos casos corresponde hacer la pregunta empírica de si acaso el fenómeno que estoy estudiando cumple con la propiedad de Markov. Por ejemplo ¿Es la evolución del clima un proceso de Markov?
#
# En segundo lugar las cadenas de Markov son un **componente fundamental de una clase de algoritmos conocidos como Markov Chain Monte Carlo** (MCMC)
#
# El objetivo de MCMC es crear sintéticamente una cadena de Markov que converge a una distribución en la cual estamos interesados y que no podemos simular de forma analítica y/o explícita
#
# MCMC es considerado una revolución en computación científica y es usado en prácticamente todos las disciplinas. 
#
#
# En esta lección estudiaremos el algoritmos de Metropolis, una de las formulaciones originales de MCMC y también uno de los [diez algoritmos más importantes del Siglo XX](https://www.andrew.cmu.edu/course/15-355/misc/Top%20Ten%20Algorithms.html)

# # Markov Chain Monte Carlo (MCMC)
#
# MCMC es una poderosa herramienta para muestrear y calcular valores esperados a partir de distribuciones complejas 
#
# En este sentido es una extensión de la idea básica de Monte Carlo que vimos en las primeras lecciones
#
# ## Monte Carlo y muestreo por importancia (IS)
#
# Sea una función $f()$ sobre una variable aleatoria con distribución $p(x)$
#
# Con Monte Carlo puedo estimar el valor esperado de esta función en base a muestras usando 
#
# $$
# \mathbb{E}[f(X)] \approx \frac{1}{N} \sum_{i=1}^N f(x_i) \quad x_i \sim p(x)
# $$
#
# Siempre y cuando yo pueda muestrear directamente de $p(x)$
#
# Si no puedo muestrear de $p(x)$ pero si es posible evaluarla, puedo recurrir a la técnica de muestreo por importancia (IS) que se define a cotninuación
#
# Sea una distribución de propuestas o distribución de importancia $q(x)$ de la cual yo puedo evaluar y además muestrear 
#
# $$
# \begin{align}
# \mathbb{E}_{x\sim p(x)}[f(X)] &= \int p(x) f(x) \,dx = \int q(x)  \frac{p(x)}{q(x)} f(x) \,dx \nonumber \\
# &= \mathbb{E}_{x\sim q(x)}\left[ \frac{p(x)}{q(x)} f(X)\right] \nonumber \\
# &\approx \frac{1}{N} \sum_{i=1}^N w_i f(x_i) \quad x_i \sim q(x) \nonumber
# \end{align}
# $$
#
# donde $w_i = \frac{p(x)}{q(x)}$ se llama la ponderación de importancia. 
#
# Una distribución de importancia correcta no sólo nos permite resolver el problema sino que tiende a tener una varianza más baja que el estimador original de Monte Carlo. No es necesario escoger una distribución de importancia que sea igual a la distribución original, pero se debe tener en cuanta que que $q(x)$ debe ser tal que $p(x)=0$ cuando $q(x)=0$ 
#

# ### Ejemplo
#
# Sea una linea de teléfono de soporte tecnológico que recibe en promedio 2 llamadas por minuto
#
# ¿Cuál es la probabilidad de que ellos tengan que esperar por lo menos 10 minutos para recibir 9 llamadas?
#
# Usemos una simulación para resolver este problema
#
# Note como el estimador basado en IS converge más rápido y con menos varianza

b = 2. # Eventos promedio por minuto
a = 9 # Cantidad de eventos
# La distribución gamma modela tiempos de espera para una cierta cantidad de eventos
p = scipy.stats.gamma(a, scale=1/b) 
# La función f en esta caso me dice 
f = lambda x: x > 10
# La función de propuesta
q = scipy.stats.norm(scale=10)
# Simulación
mc_result = []
is_result = []
true_result = 1 - p.cdf(10)
Ns = np.logspace(1, 4, num=100)
for N in Ns:
    # Monte carlo clasico
    X = p.rvs(size=int(N))
    mc_result.append(np.mean(f(X)))
    # Muestreo por importancia
    X = q.rvs(size=int(N))
    w = p.pdf(X)/q.pdf(X)
    is_result.append(np.mean(w*f(X)))
# Visualización
fig, ax = plt.subplots(figsize=(6, 3), tight_layout=True)    
ax.plot(Ns, mc_result, label='MC')
ax.plot(Ns, is_result, label='IS')
ax.axhline(true_result, c='r', ls='--', label='Real')
ax.legend()
ax.set_ylim([-0.001, true_result*3])
ax.set_xscale('log')

# ## Problemas de IS
#
# Muestreo por importancia y muestreo por rechazo me permiten calcular valores esperados de distribuciones que puedo evaluar pero no muestrear. También vimos que favorece en la disminución de la varianza
#
# Pero existen casos más complicados aun, por ejemplo 
#
# ### No puedo muestrear ni evaluar  la distribución de interés
#
# Digamos que estamos interesados en la distribución de una variable $\theta$ condicionada a un conjunto de observaciones $D$, esto corresponde al posterior $p(\theta|D)$
#
# Sólo en contadas ocasiones este posterior corresponderá a una distribución teórica como las que hemos visto en este curso
#
# Más en general tendremos
#
# $$
# p(\theta|x) = \frac{p(x|\theta) p(\theta)}{p(x)}
# $$
#
# donde $p(x|\theta)$ es la verosimilitud, $p(\theta)$ es el prior y
#
# $$
# p(x) = \int_\theta p(x, \theta) \,d\theta
# $$
#
# es la evidencia o verosimilitud marginal que no depende de $\theta$. Si la dimensionalidad de $\theta$ es grande la integral será muy difícil o derechamente imposible de calcular analiticamente. 
#
# En este caso sólo podemos evaluar la verosimilitud y el prior, es decir que podemos evaluar una distribución proporcional al posterior 
#
# $$
# p(\theta|x) \propto p(x|\theta) p(\theta)
# $$
#
# hasta la constante $1/p(x)$
#
# ### Espacios demasiado grandes
#
# Otro problema de los espacios de alta dimensionalidad es que recorrer ese espacio completo de forma independiente puede ser muy lento o de plano infactible
#

# ## ¿Cómo MCMC me salva en este caso? Intuición
#
# En MCMC en lugar de muestrear de manera iid, utilizamos una cadena de Markov que corresponde a la secuencia de pasos que damos en el espacio de alta dimensionalidad. 
#
# En la siguiente figura la distribución de interés se muestra de color rojo. En la subfigura de la izquierda usamos una distribución de importancia simple (contorno verde). Muchos de los puntos tendrán un ponderador de importancia cercano a cero. 
#
# <img src="images/is_mcmc.png" width="500">
#
# Los métodos de MCMC se basan en "diseñar" esta cadena de Markov tal que converja a la distribución complicada que nos interesa, como muestra la subfigura de la derecha
#
# Luego sólo tenemos que dejar que la cadena corra "por un tiempo largo" para que la convergencia se cumpla y finalmente usar los valores de los estados de la cadena en representación de la distribución a la cual no tenemos acceso
#
#
# Algunas definiciones:
#
# - Esta secuencia de valores se llama **traza**
# - El tiempo que demora en converger la cadena se llama **mixing time**
# - Se suele ignorar las primeras muestras de la secuencia puesto que la cadena no ha convergido. Para esto se selecciona un tiempo de **burn-in**. Luego de que se cumpla este tiempo aceptamos las muestras
#

# ## ¿Qué es diseñar una cadena de Markov?
#
# Extendiendo al caso de un estado continuo en lugar de discreto, la distribución estacionaria $\pi$ debe cumplir
#
# $$
# \int \pi(\theta_t) q(\theta_{t+1}|\theta_t) \,d\theta_t = \pi (\theta_{t+1})
# $$
#
# Diseñar la cadena de Markov consiste en encontrar las probabilidades de transición $q(\theta_{t+1}|\theta_t)$ dado que conozco $\pi$ 
#
# Notemos que esto es "al reves" de lo que hicimos en la lección pasada, que era encontrar $\pi$ dado que conozco la matriz de transición
#
# A continuación veremos veremos que no necesitamos conocer "completamente" $\pi$ para lograr esto, basta conocerlo hasta una constante

# ##  Algoritmo de Metropolis
#
# El algoritmo de Metropolis fue el primer algoritmo de tipo MCMC. Fue propuesto por <NAME>, colega de Ulam y <NAME>, [en el año 1953 para entender la transición de fase que experimetan los materiales](https://www.semanticscholar.org/paper/Equation-of-state-calculations-by-fast-computing-Metropolis-Rosenbluth/f6a13f116e270dde9d67848495f801cdb8efa25d). En el paper original sentó las bases de lo que hoy conocemos como el algoritmo de Metropolis y el algoritmo de Simulated Annealing (SA)
#
# El algoritmo de Metropolis utiliza un random walk para definir las probabilidades de transición de la cadena
#
# Sea 
#
# $$
# \theta_{t+1} = \theta_{t} + \epsilon
# $$
#
# donde $\epsilon$ se distribuye según una distribución centrada en cero y simétrica, muy tipicamente una gaussiana $\epsilon \sim \mathcal{N}(0, I\sigma_\epsilon^2)$, donde $\sigma_\epsilon$ pasa a ser un hiper parámetro del algoritmo
#
# Por definición tenemos entonces 
#
# $$
# \theta^* \sim q(\theta_{t+1}|\theta_{t}) = \mathcal{N}(\theta_{t}, I \sigma_\epsilon^2)
# $$
#
# La distribución $q$ se denomina **distribución de propuestas** y su objetivo es **proponer** un valor para  $\theta_{t+1}$ 
#
# El nuevo valor se acepta con una tasa definida como
#
# $$
# \alpha(\theta^*|\theta_{t}) = \min(1, r)
# $$
#
# donde
#
# $$
# r = \frac{ p(\theta^*)q(\theta_{t}|\theta^*) }{ p(\theta_t)q(\theta^*|\theta_{t})} = \frac{p(\theta^*)}{p(\theta_t)}
# $$
#
# donde la última equivalencia se tiene por la simetría
#
# Entonces
#
# - Si $\theta^*$ es mucho mejor que $\theta_t$ entonces se acepta
# - Si $\theta^*$ es mucho peor que $\theta_t$ entonces se rechaza
# - En caso de duda se deja al azar
#
# Respecto de $\sigma_\epsilon$
# - Si su valor es grande tendremos muchos rechazos
# - Si su valor es pequeño la difusión será lenta y podrían requerirse muchas iteraciones

# ### Formalismo
#
# El algoritmo completo es
#
#
# - Escoger una distribución de propuestas simétrica 
# - Escoger un valor inicial $\theta_0$
# - Para $n=1,2,\ldots, N$
#     - Muestrear $\theta^* \sim q(\theta_{t+1}|\theta_{t})$
#     - Muestrear $u \sim U[0, 1]$ 
#     - Luego si 
#     $$
#     u < \alpha(\theta^*|\theta_{t}) 
#     $$
#     entonces
#     $$
#     \theta_{t+1} = \theta^*
#     $$
#     de lo contrario
#     $$
#     \theta_{t+1} = \theta_{t}
#     $$
#

# ### Posteriors
#
# Notemos que si estamos interesados en un posterior, entonces
#
# $$
# r = \frac{p(\theta^*|\mathcal{D})}{p(\theta_t|\mathcal{D})} = \frac{p(\mathcal{D}|\theta^*)p(\theta^*)}{p(\mathcal{D}|\theta_t)p(\theta_t)} \frac{p(\mathcal{D})}{p(\mathcal{D})} = \frac{p(\mathcal{D}|\theta^*)p(\theta^*)}{p(\mathcal{D}|\theta_t)p(\theta_t)}
# $$
#
# Es decir que no necesitamos conocer la evidencia o verosimilitud marginal. Basta con conocer la verosimilitud y el prior

# ### Ejemplo
#
# Sea un conjunto de muestras con $N=5$
#
# $$
# \mathcal{D} = \{ 9.37, 10.18, 9.16, 11.60, 10.33 \}
# $$
#
# que corresponden a realizaciones i.i.d 
#
# $$
# X_1, X_2, \ldots, X_5|\theta \sim \mathcal{N}(\theta, \sigma^2=1)
# $$
#
# donde
#
# $$
# \theta \sim \mathcal{N}(\mu=5, \tau^2=10)
# $$
#
# y nos interesa el posterior $p(\theta|\mathcal{D})$
#
# En este caso particular el posterior si tiene una forma analítica
#
# $$
# p(\theta|\mathcal{D}) = \mathcal{N}\left ( \bar x (1- w_N) + \mu w_N , \tau_N^2 \right)
# $$
#
# donde $w_N = \tau_N^2/\tau^2$ y $\tau_N^2 = (N/\sigma^2 + 1/\tau^2)^{-1}$
#
# Intentemos simular este posterior con el algoritmo de Metropolis

# +
x = np.array([9.37, 10.18, 9.16, 11.60, 10.33])
tn2 = (len(x)/1. + 1./10)**(-1)
wn = tn2/10.

prior = lambda theta : scipy.stats.norm(loc=5, scale=np.sqrt(10)).pdf(theta)
likelihood = lambda theta : np.prod([scipy.stats.norm(loc=theta, scale=1.).pdf(x_) for x_ in x])
r = lambda ts, tt : likelihood(ts)*prior(ts)/(likelihood(tt)*prior(tt)) 

def metropolis(mix_time=5000, sigma_eps=1.):
    thetas = np.zeros(shape=(mix_time, ))
    thetas[0] = np.random.randn()
    qs = scipy.stats.norm(loc=0, scale=sigma_eps).rvs(size=mix_time)
    us = scipy.stats.uniform.rvs(size=mix_time)
    for n in range(1, mix_time):
        theta_star = thetas[n-1] + qs[n]    
        if us[n] < np.amin([1, r(theta_star, thetas[n-1])]):
            thetas[n] = theta_star
        else:
            thetas[n] = thetas[n-1]
    return thetas


# +
# %%time
burn_in = 100
thetas = metropolis(mix_time=5000, sigma_eps=1.)

fig, ax = plt.subplots(1, 2, figsize=(7, 3), tight_layout=True)
ax[0].plot(thetas)
ax[0].axhline(np.mean(x)*(1-wn) + 5*wn, c='r', ls='--', lw=2, alpha=0.5)
ax[1].hist(thetas[burn_in:], density=True, bins=10)
t_plot = np.linspace(np.amin(thetas[burn_in:]), 
                     np.amax(thetas[burn_in:]), num=100)
ax[1].plot(t_plot, scipy.stats.norm(loc=np.mean(x)*(1-wn)+5*wn,
                                    scale=np.sqrt(tn2)).pdf(t_plot), 
           c='r', lw=2, ls='--', alpha=0.5);
# -

# **Propuesto**
#
# - Estudie como cambian los resultados con $\sigma_\epsilon \in \{0.01, 1, 100\}$
# - Estudie como cambian los resultados con distintos valores de $\theta_0$

# ## Algoritmo de Metropolis-Hastings
#
# El algoritmo de Metropolis Hastings es una generalización del algoritmo de Metropolis para el caso donde la distribución de propuestas ya no es simétrica por lo que
#
# $$
# r = \frac{ p(\theta^*)q(\theta_{t}|\theta^*) }{ p(\theta_t)q(\theta^*|\theta_{t})}
# $$
#
# El algoritmo procede de forma idéntica al caso anterior



StarcoderdataJupyterScriptsDedupFiltered