09-max-verosimilitud.Rmd

# Estimación por máxima verosimilitud {#S:max-verosimilitud}

```{r setup, include=FALSE, message=FALSE}
library(tidyverse)
library(patchwork)
library(broom)
source("R/funciones_auxiliares.R")
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning=FALSE, fig.align = 'center', fig.width = 5, fig.height=3)
comma <- function(x) format(x, digits = 2, big.mark = ",")
theme_set(theme_minimal())
```


Los ejemplos que hemos visto han sido todos de estimadores *plug-in* (o por sustitución): 
si queremos saber una cantidad poblacional, y tenemos una muestra dada, entonces
calculamos la estadística de interés *como si la muestra fuera la población*. Por ejemplo,
para estimar la mediana poblacional usamos la mediana muestral, si queremos estimar
la media poblacional usamos la media muestral, y así sucesivamente. Estos estimadores
usualmente dan resultados razonables (pero hay que checar usando muestra bootstraps, por ejemplo,
y pensar lo que estamos haciendo).

Cuando sabemos más acerca de la población y usamos un modelo teórico es posible
hacer más: dependiendo de qué cantidades se quieren estimar, podemos construir estimadores
que sean *óptimos* en algún sentido siempre y cuando se cumplan los supuestos teóricos, como veremos ahora. 

Por ejemplo: ¿deberíamos estimar el centro de una distribución simétrica con la
media o con la mediana, o quizá con una media recortada?

En esta parte construiremos la teoría básica de estimación cuando trabajamos con 
modelos teóricos conocidos. El objetivo es entender **las ideas básicas** de estos procedimientos, 
y cómo evaluar sus resultados.

**Recordatorio**: las ventajas de usar modelos teóricos para describir distribuciones de
datos está en que es posible comprimir más eficientemente la información, es posible construir
modelos más complejos juntando varios de estos modelos y de sus dependencias, y de que es
posible hacer más teoría útil que nos guíe. La desventaja es que es necesario que esos supuestos teóricos sean razonables.


## Introducción a estimación por máxima verosimilitud {-}

Uno de los procedimientos más estándar en esta situación es el **método 
de máxima verosimilitud**. Los estimadores de máxima verosimilitud tienen propiedades
convenientes, y dan en general resultados razonables siempre y cuando los supuestos
sean razonables. 

**Máxima verosimilitud es un proceso intuitivo, y consiste en aprender o estimar valores de parámetros desconocidos suponiendo para los datos su explicación más probable**. Para
esto, usando supuestos y modelos requeriremos calcular la probabilidad de un conjunto
de observaciones.

**Ejemplo.** Adaptado de [@Chihara]. Supongamos que una máquina produce dos tipos de bolsas de 25 galletas: la mitad de las veces produce una bolsa con 5 galletas de avena y 20 de chispas de chocolate, y
la otra mitad produce galletas con 23 galletas de avena y 2 de chispas de chocolate. 

Tomamos una bolsa, y no sabemos qué tipo de bolsa es (parámetro desconocido). Extraemos al azar una de las galletas, y es de chispas de chocolate (observación).

Por máxima verosimilitud, inferimos que la bolsa que estamos considerando tiene 5 galletas de avena. Esto es porque es más probable observar una galleta de chispas en las bolsas que contienen 5 galletas de avena que en las bolsas que contienen 23 galletas de avena. Podemos cuantificar la probabilidad que "acertemos" en nuestra inferencia.



Cómo se aprecia en el ejemplo anterior, el esquema general es:

1. Existe un proceso del que podemos obtener observaciones de algún sistema
o población real.
2. Tenemos un modelo probabilístico que dice cómo se producen esas observaciones a partir del sistema o población real.
3. Usualmente este modelo tiene algunas cantidades que no conocemos, que rigen el proceso y cómo se relaciona el proceso con las observaciones.

Nuestro propósito es:

4. Extraemos observaciones del proceso  
$$x_1, x_2, \ldots, x_n.$$
5. Queremos aprender de los parámetros desconocidos del proceso para calcular cantidades de interés acerca del sistema o población real


En principio, los modelos que consideramos pueden ser complicados y tener varias partes o parámetros. Veamos primero un ejemplo clásico con un solo parámetro, y cómo lo resolveríamos
usando máxima verosimilitud.

**Nota**: Cuando decimos *muestra* en general nos referimos a observaciones
independientes obtenidas del mismo proceso (ver la sección \@ref(S:distribucion-muestreo) para ver qué 
significa que sea independientes). Este esquema es un supuesto
que simplifica mucho los cálculos, como discutimos antes. Muchas veces este supuesto
sale del diseño de la muestra o del estudio, pero en todo caso es importante considerar
si es razonable o no para nuestro problema particular.

```{block, type = 'mathblock'}
Denotemos por $f(x; \theta)$ la función de densidad para una variable aleatoria continua con párametro 
asociado $\theta.$ Denotamos por $X_1, \ldots, X_n,$ una muestra aleatoria de $n$ observaciones de esta 
distribución y por $x_1, \ldots, x_n$ los valores observados de esta muestra aleatoria.
```


**Ejemplo.** Supongamos que queremos saber qué proporción de registros de una base de datos tiene algún error menor de captura. No podemos revisar todos los registros, así que tomamos una muestra
de 8 registros, escogiendo uno por uno al azar de manera independiente. Revisamos los
8 registros, y obtenemos los siguientes datos:

$$x_1 = 0, x_2 = 1, x_3 = 0, x_4 = 0, x_5 =1, x_6 =0, x_7 =0, x_8 =0$$

donde 1 indica un error menor. Encontramos dos errores menores. ¿Cómo estimamos el número de
registros con errores leves en la base de datos?

Ya sabemos una respuesta razonable para nuestro estimador puntual, que sería $\hat{p}=2/8=0.25$. 
Veamos cómo se obtendría por máxima verosimilitud.

Según el proceso con el que se construyó la muestra, debemos dar una probabilidad de observar
los 2 errores en 8 registros. Supongamos que en realidad existe una proporción $p$
de que un registro tenga un error. Entonces calculamos

Probabilidad de observar la muestra:

$$P(X_1 = 0, X_2 = 1, X_3 = 0, X_4 = 0, X_5 =1, X_6 =0, X_7 =0, X_8 =0)$$

es igual a

$$P(X_1 = 0)P(X_2 = 1)P(X_3 = 0)P( X_4 = 0)P(X_5 =1)P(X_6 =0)P(X_7 =0)P(X_8 =0)$$
pues la probabilidad de que cada observación sea 0 o 1 no depende de las observaciones restantes (la muestra se extrajo de manera independiente).

Esta última cantidad tiene un parámetro que no conocemos: la proporcion $p$ de registros
con errores. Así que lo denotamos como una cantidad desconocida $p$. Nótese
entonces que $P(X_2=1) = p$, $P(X_3=0) = 1-p$ y así sucesivamente, así que la cantidad
de arriba es igual a

$$(1-p)p(1-p)(1-p)p(1-p)(1-p)(1-p) $$

que se simplifica a 

$$ \mathcal{L}(p) = p^2(1-p)^6$$

Ahora la idea es **encontrar la p que maximiza la probabilidad de lo que observamos**. En
este caso se puede hacer con cálculo, pero vamos a ver una gráfica de esta función y
cómo resolverla de manera numérica.

```{r}
verosimilitud <- function(p){
  p^2 * (1-p)^6
}
dat_verosim <- tibble(x = seq(0,1, 0.001)) %>% mutate(prob = map_dbl(x, verosimilitud))
ggplot(dat_verosim, aes(x = x, y = prob)) + geom_line() +
  geom_vline(xintercept = 0.25, color = "red") +
  xlab("p")
```


Nótese que esta gráfica:

- Depende de los datos, que pensamos fijos.
- Cuando cambiamos la $p$, la probabilidad de observar la muestra cambia. Nos interesa
ver las regiones donde la probabilidad es relativamente alta.
- El máximo está en 0.25.
- Así que el estimador de máxima verosimilitud es $\hat{p} = 0.25$, **que es también
el estimador usual de plugin** en este caso.


Para uniformizar la notación con el caso continuo que veremos más adelante, usaremos
la notación 

$$P(X=x) = f(x)$$
donde $f$ es la función de densidad (en este caso, función de masa de probabilidad) de $X$. Si esta función depende
de un parámetro, escribimos $$f(x ;\theta)$$


```{block verosimilitud, type = 'mathblock'}
**Definición.** Sean $X_1, \ldots, X_n$ una muestra de una densidad  $f(x; \theta)$
y sean $x_1,x_2,\ldots, x_n$ los valores observados. 

La *función de verosimilitud* del párametro de interés $\theta$ está definida por 
\begin{align}
\mathcal{L}(\theta; x_1, \ldots, x_n) = \prod_{i = 1}^n f(x_i; \theta). 
\end{align}

Esta función nos dice qué tan creible es el valor del parámetro $\theta$ dada la muestra observada. 
A veces también la denotamos por $\mathcal{L}_n(\theta)$. 
```


Ahora definimos qué es un estimador de máxima verosimilitud.


```{block mle, type = 'mathblock'}
**Definición.** Un estimador de máxima verosimilitud lo denotamos por $\hat \theta_{\textsf{MLE}}$ y es un valor que satisface 
\begin{align}
\hat \theta_{\textsf{MLE}} =  \underset{\theta \, \in \, \Theta}{\arg\max}\,  \mathcal{L}(\theta; x_1, \ldots, x_n),
\end{align}
donde $\theta$ denota el espacio parametral. Es decir, el espacio válido de búsqueda congruente con la definición del modelo. 
```

```{block, type = 'ejercicio'}

- Considera el caso de una normal con media y varianza desconocidas. ¿Cuáles son los espacios parametrales para efectuar $\mathsf{MLE}$?
- Considera el caso de una Binomial con parámetros $N$ y $p$ desconocidos. ¿Cuáles son los espacios parametrales para la búsqueda del $\mathsf{MLE}$?
```


Obsérvese que para construir la verosimilitud y en consecuencia buscar por
estimadores de máxima verosimlitud necesitamos:

- Un modelo teórico de cómo es la población con parámetros e  
- Información de cómo se extrajo la muestra,  

y entonces podemos resolver nuestro problema de estimación 
convirtiéndolo en uno de optimización.

Probamos esta idea con un proceso más complejo.


**Ejemplo.** Supongamos que una máquina puede estar funcionando correctamente o no en cada
corrida. Cada corrida se producen 500 productos, y se muestrean 10 para
detectar defectos. Cuando la máquina funciona correctamente, la tasa de defectos
es de 3%. Cuando la máquina no está funcionando correctamente la tasa de defectos
es de 20%

Supongamos que escogemos al azar 11 corridas, y obervamos los siguientes
número de defectuosos:

$$1, 0, 0, 3 ,0, 0, 0, 2, 1, 0, 0$$

La pregunta es: ¿qué porcentaje del tiempo la máquina está funcionando correctamente?

Primero pensemos en una corrida. La probabilidad de observar una sucesión particular de 
$r$ defectos es

$$0.03^r(0.97)^{(10-r)}$$
cuando la máquina está funcionando correctamente.

Si la máquina está fallando, la misma probabilidad es 

$$0.2^r(0.8)^{(10-r)}.$$

Ahora supongamos que la máquina trabaja
correctamente en una proporción $p$ de las corridas. Entonces la probabilidad
de observar $r$ fallas se calcula promediando (probabilidad
total) sobre las probabilidades de que la máquina esté funcionando
bien o no:

$$0.03^r(0.97)^{(10-r)}p + 0.2^r(0.8)^{(10-r)}(1-p)$$
Y esta es nuestra función de verosimilitud para una observación.

Suponemos que las $r_1,r_2, \ldots, r_{11}$ observaciones son independientes
(por ejemplo, después de cada corrida la máquina se prepara de una manera estándar,
y es como si el proceso comenzara otra vez). Entonces tenemos que multiplicar
estas probabilidades para cada observación $r_1$:


```{r}
calc_verosim <- function(r){
  q_func <- 0.03^r*(0.97)^(10-r)
  q_falla <- 0.2^r*(0.8)^(10-r)
  function(p){
    #nota: esta no es la mejor manera de calcularlo, hay 
    # que usar logaritmos.
    prod(p*q_func + (1-p)*q_falla)
  }
}
verosim <- calc_verosim(c(1, 0, 0, 3, 0, 0, 0, 2, 1, 0, 0))
verosim(0.1)
```

```{r}
dat_verosim <- tibble(x = seq(0,1, 0.001)) %>% mutate(prob = map_dbl(x, verosim))
ggplot(dat_verosim, aes(x = x, y = prob)) + geom_line() +
  geom_vline(xintercept = 0.8, color = "red") +
  xlab("prop funcionado")
```

Y nuestra estimación puntual sería de alrededor de 80%.


## Máxima verosimilitud para observaciones continuas {-}

Cuando las observaciones $x_1,\ldots, x_n$ provienen de una distribución
continua, no tiene sentido considerar $P(X = x_i)$, pues siempre es igual
a cero. 

Sin embargo, podemos
escribir para pequeños valores $\epsilon \ll 1$
\begin{align}
  P(x - \epsilon < X < x + \epsilon | \theta) = \int_{x - \epsilon}^{x + \epsilon} f(t; \theta) \, \text{d} t \approx 2 \epsilon f(x; \theta),
\end{align}
donde $f(x; \theta)$ es la función de densidad de $X.$ Por lo tanto, 
\begin{align}
\begin{split}
  P&(x_1 - \epsilon < X_1 < x_1 + \epsilon, \ldots, x_n - \epsilon < X_n < x_n + \epsilon | \theta) \\
  &= \prod_{i = 1}^n P(x_i - \epsilon < X_i < x_i + \epsilon | \theta) \\
  &= \prod_{i = 1}^n 2 \epsilon f(x_i; \theta) = (2\epsilon)^n \prod_{i = 1}^n f(x_i; \theta).
\end{split}
\end{align}

Notemos que si $\epsilon \rightarrow 0$ la ecuación rápidamente converge a cero. Pero para pequeños 
valores de $\epsilon$ la ecuación que nos interesa es proporcional a $\prod_{i = 1}^n f(x_i; \theta).$

De esta forma, nuestra definición de máxima verosimilitud y estimadores
de máxima verosimilitud es la misma para el caso continuo (verifica las definiciones
de la sección anterior).

**Ejemplo.** Supongamos que tenemos una muestra $x_1\ldots, x_n$ extraidas de una distribución
exponencial con tasa $\lambda>0$ donde no conocemos $\lambda$. ¿Cuál es 
el estimador de máxima verosimilitud de $\lambda$?

Para $\lambda>0$, tenemos que

$${\mathcal L}(\lambda) = \prod_{i=1}^n \lambda e^{\lambda x_i}$$
de modo que

$${\mathcal L}(\lambda) = \lambda^n e^{-\lambda \sum_{i=1}^nx_i} = \lambda^n e^{-n\lambda\bar{x}} = e^{n(\log\lambda - \lambda\bar{x})}$$
Que podemos maximizar usando cálculo para obtener 
$\hat{\lambda}_{\mathsf{ML}} = \frac{1}{\bar{x}}$ (demuéstralo). Discute
por qué esto es intuitivamente razonable: ¿cuál es
el valor esperado de una exponencial con parámetro $\lambda$?


## Aspectos numéricos {-}

Encontrar el estimador de máxima verosimilitud ($\textsf{MLE}$) es automático en la mayoría de los casos.
En teoría, podemos reutilizar la misma rutina numérica para encontrar el estimador sin ninguna ayuda de la analista. 
Esto contrasta con otras técnicas de estimación en donde se requieren cálculos y manipulación de ecuaciones. 

Sin embargo, hay situaciones que se pueden evitar de manera general. Por ejemplo, cuando calculamos la verosimilitud arriba, nótese que estamos multiplicando
números que pueden ser muy chicos (por ejemplo $p^6$, etc). Esto puede producir
desbordes numéricos fácilmente. Por ejemplo para un tamaño de muestra de 1000, podríamos
tener que calcular

```{r}
p <- 0.1
proba <- (p ^ 800)*(1-p)^200
proba
```

En estos casos, es mejor hacer los cálculos en escala logarítmica. El logaritmo
convierte productos en sumas, y baja exponentes multiplicando. Si calculamos en escala
logaritmica la cantidad de arriba, no tenemos problema:

```{r}
log_proba <- 800 * log(p) + 200 * log(1-p)
log_proba
```

Ahora notemos que 

- Maximizar la verosimilitud **es lo mismo** que maximizar la log-verosimilitud,  pues el logaritmo es una función creciente. Si $x_{\max}$ es el máximo de $f$, tenemos que $f(x_{\max})>f(x)$ para cualquier $x$, entonces tomando logaritmo,
$$\log(f(x_{max}))>\log(f(x)),$$ para cualquier $x.$ Pues el logaritmo respeta la desigualdad por ser creciente.
- Usualmente usamos la log-verosimilitud para encontrar el estimador de máxima verosimilitud.
- Hay razónes teóricas y de interpretación por las que también es conveniente hacer esto.

```{block log-like, type = 'mathblock'}
**Definición.** La log-verosimilitud la denotamos usualmente por 
$$\ell_n(\theta) = \log \left(\mathcal{L}_n(\theta)\right),$$
donde hemos suprimido la dependencia en la muestra por conveniencia.
```

```{block, type = 'ejercicio'}
Considera una muestra de variables aleatorias Gaussianas con media $\mu$ y varianza $\sigma^2.$ Escribe la 
verosimilitud para una muestra de tamaño $n,$ y después escribe la función de log-verosimilitud. ¿Que interpetación 
le das a la ecuación resultante? ¿La has visto en algunos otros ejemplos en secciones anteriores?  
*Pista.* Recuerda la sección de regresión local.  
```



**Ejemplo.** En nuestro primer ejemplo,

```{r}
log_verosimilitud <- function(p){
  2*log(p) + 6*log(1-p)
}
dat_verosim <- tibble(x = seq(0,1, 0.01)) %>% mutate(log_prob = map_dbl(x, log_verosimilitud))
ggplot(dat_verosim, aes(x = x, y = log_prob)) + geom_line() +
  geom_vline(xintercept = 0.25, color = "red") +
  xlab("p")
```

Obtenemos el mismo máximo. Podemos incluso resolver numéricamente:

```{r, warning = FALSE}
solucion <- optim(p = 0.5, log_verosimilitud, control = list(fnscale = -1))
solucion$par
```

Y en nuestro segundo ejemplo:

```{r}
calc_log_verosim <- function(r){
  q_func <- 0.03^r*(0.97)^(10-r)
  q_falla <- 0.2^r*(0.8)^(10-r)
  function(p){
    #nota: esta no es la mejor manera de calcularlo, hay 
    # que usar logaritmos.
    sum(log(p*q_func + (1-p)*q_falla))
  }
}
log_verosim <- calc_log_verosim(c(1, 0, 0, 3, 0, 0, 0, 2, 1, 0, 0))
log_verosim(0.1)
```

```{r}
dat_verosim <- tibble(x = seq(0,1, 0.001)) %>% mutate(log_verosimilitud = map_dbl(x, log_verosim))
ggplot(dat_verosim, aes(x = x, y = log_verosimilitud)) + geom_line() +
  geom_vline(xintercept = 0.775, color = "red") +
  xlab("prop funcionado")
```



En la función de verosimilitud los datos están fijos y optimizamos con respecto a
las variables que no conocemos. ¿Cómo construimos esta función en general, para cualquier
tamaño de muestra y cualquier número de registros detectados como erróneos?

Nótese que la verosimilitud la consideramos **función de los parámetros**,
donde **los datos están fijos**.

Podemos construir una función que genera la función de verosimilitud dependiendo de los datos. En nuestro primer ejemplo de muestras de registros erróneos,
podríamos construir una función que genera la log verosimilitud dependiendo
del tamaño de muestra y del número de errores encontrado:

```{r}
construir_log_verosim <- function(n, n_err){
  # n es tamaño de muestra 
  # n_err el número de errores detectados (datos)
  n_corr <- n - n_err
  log_verosim <- function(p){
    n_err * log(p) + n_corr * log(1-p)
  }
}
```

Cuando fijamos $n$ y $n_{\textsf{err}}$, esta función genera otra función, la log verosimilitud, que es la que queremos optimizar.

Supongamos entonces que sacamos 20 registros al azar y observamos 10 incorrectos. La función
de verosimilitud es

```{r}
log_vero <- construir_log_verosim(20, 10)
```


```{r}
tibble(x = seq(0,1,0.001)) %>% 
  mutate(log_ver = log_vero(x)) %>% 
  ggplot(aes(x = x, y = log_ver)) + 
    geom_line() +
    geom_vline(xintercept = 0.5, color = 'red')
```

**Ejemplo.** Supongamos que en una población de transacciones hay un porcentaje $p$ (desconocido) 
que son fraudulentas. Tenemos un sistema de clasificación humana que que marca transacciones como sospechosas. 
Con este sistema hemos medido que la proporción de transacciones normales que son marcadas como sospechosas es de 0.1%, y que la proporción de transacciones fraudulentas que son marcadas
como sospechosas es de 98%. Supongamos que extraemos una muestra de 2000 transacciones, de manera que todas ellas tiene la misma probabilidad de ser fraudulentas. El sistema de clasificación marca 6 transacciones como fraudulentas. ¿Cómo estimamos la proporción de transacciones fraudulentas en la población?

Solución: sea $p$ la proporción de transacciones fraudulentas. Entonces la probabilidad
de que una transacción sea marcada como sospechosa es (proba total):

$$0.98p + 0.001(1-p)$$

Pues tenemos que contar 98% de la proporción $p$ de fraudulentas 
(correctamente detectadas) más 0.1% de la proporción $(1-p)$ de fraudulentas.
Escribimos entonces nuestra función de verosimilitud

```{r}
crear_log_verosim <- function(n, n_sosp){
  # devolver la función log verosimilitud 
  log_verosimilitud_pct <- function(pct){
    # sup que pct es la proporcentaje de fraudes,
    # que es el parámetro que queremos estimar
    prob_sosp <- 0.98 * pct / 100 + 0.001 * (1 - pct / 100)
    log_prob <- n_sosp * log(prob_sosp) + (n - n_sosp) * log(1- prob_sosp)
    log_prob
  }
  log_verosimilitud_pct
}
```

La verosimilitud es una función de $p$.

```{r}
log_verosim <- crear_log_verosim(n = 2000, n_sosp = 4)
```

A continuación la mostramos de manera gráfica.

```{r, echo = FALSE}
tibble(x = seq(0,10,0.001)) %>% 
  mutate(log_ver = log_verosim(x)) %>% 
  ggplot(aes(x = x, y = log_ver)) + 
    geom_line() +
    xlab("Porcentaje fraudulentas")
```

No se ve muy claro dónde ocurre el máximo, pero podemos ampliar cerca de cero la 
misma gráfica:

```{r, echo = FALSE}
tibble(x = seq(0,0.5,0.001)) %>% 
  mutate(log_ver = log_verosim(x)) %>% 
  ggplot(aes(x = x, y = log_ver)) + 
    geom_line() +
    xlab("Porcentaje fraudulentas")
```



- Vemos que alrededor de 0.1% maximiza la probabilidad de haber observado 4 
transacciones sospechosas. 
- Notamos sin embargo que varios valores alrededor de este valor tienen probabilidad similar,
así que también son consistentes con los datos (por ejemplo, valores como 0.05 o 0.15 tienen probabilidad similar). Tendremos que considerar esto para evaluar la incertidumbre en nuestra estimación.
- Obsérvese adicionalmente que si no tomáramos en cuenta las probabilidades de falsos
negativos y falsos positivos la estimación simple daría $6/2000 = 0.003$ (0.3%), que es
tres veces más grande que nuestra estimación puntual por máxima verosimilitud.


**Ejemplo.** Este es un ejemplo donde mostramos que cuando el soporte
de las densidades teóricas es acotado hay que tener cuidado en la definición de la 
verosimilitud. En particular, el soporte de la variable aleatoria es el párametro de interés. Supongamos
que nuestros datos son generados por medio de una distribución uniforme en el intervalo $[0,b].$
Contamos con una muestra de $n$ observaciones generadas 
de manera independiente $X_i \sim U[0,b]$ para $i= 1, \ldots, n.$
Sin embargo, no conocemos el valor de $b$. 

¿Cómo es la función de log verosimilitud ${\mathcal L}_n(b)$ para este caso? Nótese
que cuando el parámetro $b$ es menor que alguna $x_i$, tenemos que
${\mathcal L}_n(b) = 0$: la verosimilitud es cero si tomamos una $b$ más chica
que algún dato, pues este valor es incosistente del todo con los datos observados.
En otro caso,

$${\mathcal L}_n(b) = \frac{1}{b^n},$$
pues la función de densidad de una uniforme en $[0,b]$ es igual a $1/b$ en 
el intervalo $[0,b]$, y 0 en otro caso. Podemos escribir entonces:


```{r}
crear_verosim <- function(x){
  n <- length(x)
  verosim <- function(b){
    indicadora <- ifelse(all(x <= b), 1, 0)
    indicadora / b^n
  }
}
```



Ahora podemos hacer máxima verosimilitud para un ejemplo:

```{r}
set.seed(234)
x <- runif(10, 0, 3)
verosim <- crear_verosim(x)
res_opt <- optimize(verosim, c(-1000, 1000), maximum = TRUE)
res_opt$maximum
```
Y nótese que, como esperaríamos, este valor es el máximo de la muestra:
```{r}
max(x)
```

La gráfica de la función de verosimilitud es:

```{r}
tibble(b = seq(-1, 5, 0.001)) %>% 
  mutate(verosim_1 = map_dbl(b, ~ verosim(.x))) %>% 
ggplot() + 
  geom_line(aes(x = b, y = verosim_1)) +
  geom_rug(data = tibble(x = x), aes(x = x), colour = "red")
```

Podemos escribir en una fórmula como:

\begin{align}
\mathcal{L}(b; x_1, \ldots, x_n) = \prod_{i = 1}^n 1_{[0,b]}(x_i) \frac1b. 
\end{align}

Y podríamos resolver analíticamente como sigue:

Si consideramos 
$$ \hat b_{\textsf{MLE}} = x_{\max} = \max\{x_i\},$$
notemos que cualquier valor observado necesariamente satisface 
$$x_i \leq \hat b_{\textsf{MLE}},$$ 
y por lo tanto todas las funciones indicadoras están *encendidas*. El valor de la verosimilitud es igual a 
$$\mathcal{L}(\hat b_{\textsf{MLE}}; x_1, \ldots, x_n) = \left(\frac{1}{x_{\max}}\right)^n \geq \left (\frac1b\right )^n$$ 
para cualquier $b\geq x_{\max}$. Como la verosimilitud para $b<x_{\max}$ es igual
a cero, esto demuestra que el máximo de la muestra es el estimador de máxima
verosimilitud de $b$.

**Observación.** Este ejemplo también tiene dificultades numéricas, pues
la verosimilitud presenta discontinuidades y regiones con derivada igual a cero, y
la mayoria de los algoritmos numéricos no tienen garantías buenas de 
covergencia al máximo en estos casos.
Si aplicamos sin cuidado descenso en gradiente, por ejemplo, podríamos comenzar
incorrectamente en un valor $b_0 < x_{\max}$ y el algoritmo no avanzaría
al máximo. 



### El método de momentos {-}

Un método alternativo para estimación de parámetros es el método de momentos ($\textsf{MOM}$). Para esto 
creamos un sistema de ecuaciones de cardinalidad igual al número de parámetros a estimar. 
Es decir, consideramos $\theta \in \mathbb{R}^p;$ los momentos teoricos
\begin{align}
  m_k(\theta) = \mathbb{E}_f[X^k] = \int_{\mathcal{X}} x^k \, f(x; \theta) \, \text{d} x;
\end{align}
los momentos empíricos 
\begin{align}
  \hat m_k(\theta) = \frac1n \sum_{i = 1}^n X_i^k;
\end{align}
y creamos el siguiente sistema de ecuaciones: 
\begin{align}
  m_k(\theta) = \hat m_k(\theta), \qquad k = 1, \ldots, p.
\end{align}
Este sistema explota la aproximación $m_k(\theta) \approx \hat m_k(\theta),$ cuya justificación está determinada
por la *Ley de los grandes números.*

**Ejemplo.** *(continuación).* Para el caso de  $n$ observaciones $X_i \sim U[0,b],$ el método de momentos 
arroja un estimador de la forma
$$\hat b_{\textsf{MOM}} = 2 \bar X_n.$$
```{block , type = 'ejercicio'}
Considera el caso de una muestra de tamaño $n = 3,$ con observaciones $X_1 = X_2 = 1$ y $X_3= 7$. ¿Cuál es el estimador $\textsf{MOM}$? ¿Qué implicaciones tiene sobre la observación $X_3$?
```


## Máxima verosimilitud para más de un parámetro {-}

Si nuestro modelo contiene más de un parámetro desconocido podemos también usar
máxima verosimilitud. En este caso, optimizamos sobre todos los parámetros usando
cálculo o alguna rutina numérica

**Ejemplo.** Supongamos que en una población de estudiantes tenemos dos tipos: unos llenaron un
examen de opción múltiple al azar (1 de 5), y otros contestaron las preguntas intentando
sacar una buena calificación. Suponemos que una vez que conocemos el tipo de 
estudiante, todas las preguntas tienen la misma probabilidad de ser contestadas
correctamente, de manera independiente. El modelo
teórico está representado por la siguiente simulación:

```{r}
sim_formas <- function(p_azar, p_corr){
  tipo <- rbinom(1, 1, 1 - p_azar)
  if(tipo==0){
    # al azar
    x <- rbinom(1, 10, 1/5)
  } else {
    # no al azar
    x <- rbinom(1, 10, p_corr)
  }
  x
}
```

Y una muestra se ve como sigue:

```{r}
set.seed(12)
muestra <- map_dbl(1:200, ~ sim_formas(0.3, 0.75))
qplot(muestra)
```

Supongamos que no conocemos la probabildad de contestar correctamente  ni la
proporción de estudiantes que contestó al azar. ¿Como estimamos estas dos cantidades?

Escribimos la verosimilitud:

```{r}
crear_log_p <- function(x){

  log_p <- function(pars){
    p_azar = pars[1]
    p_corr = pars[2]
    sum(log(p_azar * dbinom(x, 10, 1/5) + (1 - p_azar) * dbinom(x, 10, p_corr)))
  }  
  log_p
}
```

Creamos la función de verosimilitud con los datos

```{r}
log_p <- crear_log_p(muestra)
```

y optimizamos

```{r}
res <- optim(c(0.5, 0.5), log_p, control = list(fnscale = -1))
res$par
```

En este caso, obtenemos estimaciones razonables de ambos parámetros. 
Nota: dependiendo de los datos, este problema
puede estar mal condicionado. Por ejemplo, ¿qué pasa si la probabilidad de acertar
cuando se contesta bien está cercano al azar?

**Ejemplo.** (Tomado de [@zuev]). Considera el caso de $n$ muestras iid de un modelo Gaussiano. Es decir, 
$X_1, \ldots, X_n \sim \mathsf{N}(\mu, \sigma^2).$ Consideremos que ambos parámetros son desconocidos y nos gustaria 
encontrar el $\textsf{MLE}$. Para este problema denotamos $\theta \in \mathbb{R}^2$, donde $\theta_1 = \mu$ y $\theta_2 = \sigma^2.$ 

La función de verosimiltud se puede calcular (ignorando algunas constantes multiplicativas) como 
\begin{align}
\mathcal{L}_n(\theta) &= \prod_{i = 1}^n \frac{1}{\sigma} \, \exp\left( - \frac{(x_i - \mu)^2}{2\sigma^2}\right) \\
  &= \theta_2^{-\frac{n}{2}}\exp\left( - \frac{1}{2 \theta_2} \sum_{i = 1}^n (x_i - \theta_1)^2 \right).
\end{align}

A continuación mostramos la representación gráfica de la función de verosimilitud de este ejemplo. 
Notamos lo mismo que para los ejemplos anteriores. Conforme más datos tenemos, más nos acercamos a los valores 
reales que no conocemos. 

```{r, echo=FALSE}
knitr::include_graphics('images/mle-normal.png', auto_pdf = TRUE )
```

En las siguientes secciones veremos los métodos clásicos para encontrar el $\textsf{MLE}$. Y mas adelante, estudiaremos los propiedades teóricas que hacen a estos estimadores tan útiles y prácticos en aplicaciones. 

**Ejemplo**. Como ejercicio, podemos encontrar los estimadores de máxima
verosimilitud cuando tenemos una muestra $X_1, \ldots, X_n \sim \mathsf{N}(\mu, \sigma^2).$
(puedes derivar e igualar el cero para encontrar el mínimo). También podemos resolver numéricamente,
por ejemplo:

Supongamos que tenemos la siguiente muestra:

```{r}
set.seed(41852)
muestra <- rnorm(150, mean = 1, sd = 2)
```


La función generadora de la log verosimilitud para una muestra es (ve la expresión
del ejercicio anterior y calcula su logaritmo), y generamos la función de verosimilitud
para nuestra muestra:

```{r}
crear_log_p <- function(x){
  log_p <- function(pars){
    media = pars[1]
    desv_est = pars[2]
    n <- length(x)
    # ve la ecuación del ejercicio anterior
    z <- (x - media) / desv_est
    log_verosim <- -(log(desv_est) +  0.5 * mean(z^2))
    log_verosim
  }  
  log_p
}
log_p <- crear_log_p(muestra)
```


Ahora optimizamos (checa que el método converge):

```{r}
res <- optim(c(0, 0.5), log_p, control = list(fnscale = -1, maxit = 1000), method = "Nelder-Mead")
res$convergence
est_mv <- tibble(parametro = c("media", "sigma"), estimador = res$par) %>% 
  column_to_rownames(var = "parametro")
est_mv
```

Verifica que el estimador de la media y de la desviación estándar
es el que esperábamos (y que puedes derivar analíticamente):

```{r}
n <- length(muestra)
sd_n <- function(x) sqrt( mean((x - mean(x))^2))
c(media = mean(muestra), sigma = sd_n(muestra)) %>% round(4)
```

La siguiente pregunta qué nos interesa hacer es: ¿cómo estimamos la variabilidad
de estos estimadores? Más adelante veremos una respuesta basada en teoría, pero
también podemos resolver este problema usando el bootstrap.