---
title: "Entregable 2 - Estadística"
subtitle: "Grau de Matemàtiques, Curs 2020-21"
author: "Vilanova Martínez, Adrià"
output: html_document
---

```{r message=FALSE}
# Libraries
library(EnvStats)
```

## Configuració
Aquí es pot configurar la pràctica depenent del que es pregunti:

```{r}
# Si |generate_random| és TRUE, en comptes d'utilitzar la mostra donada a la
# variable |input|, es genera una mostra aleatòria seguint la distribució
# binomial negativa de mida |random_n|, amb paràmetres |random_size| i
# |random_prob|
generate_random = FALSE

# Dades de la mostra (s'usen si |generate_random| és FALSE)
input = data.frame(
  value = seq(0, 29),
  freq = c(0, 0, 2, 5, 11, 34, 68, 77, 139, 152, 192, 211, 190, 170, 170, 155, 109, 98, 69, 59, 31, 23, 20, 6, 3, 1, 1, 2, 2, 0) # Adri
  #freq = c(4, 26, 60, 85, 155, 177, 186, 197, 213, 202, 148, 123, 116, 74, 58, 55, 40, 31, 11, 11, 5, 11, 2, 6, 2, 1, 1, 0, 0, 0) # Pedro
  #freq = c(0, 5, 20, 31, 71, 97, 143, 155, 202, 174, 193, 167, 166, 141, 115, 84, 59, 54, 37, 28, 20, 15, 8, 6, 5, 1, 2, 0, 0, 1) # Paula
)

# Paràmetres de la mostra aleatòria (s'usen si |generate_random| és TRUE)
random_size = 10
random_prob = 0.3
random_n = 2000

# == Adri ==
# Probabilitat de la regió de predicció de la pregunta 11:
prob_11 = 0.99

# Paràmetres a, b i r de l'apartat c:
a_c = -0.3
b_c = -0.5
r_c = 43

# == Pedro ==
# Probabilitat de la regió de predicció de la pregunta 11:
#prob_11 = 0.95

# Paràmetres a, b i r de l'apartat c:
#a_c = 0.3
#b_c = 2
#r_c = 7

# == Paula ==
# Probabilitat de la regió de predicció de la pregunta 11:
#prob_11 = 0.9

# Paràmetres a, b i r de l'apartat c:
#a_c = 0.2
#b_c = 0.8
#r_c = 15
```

```{r}
if (generate_random == FALSE) {
  dd = data.frame(
    v = rep(input$value, input$freq)
  )
} else {
  dd = data.frame(
    v = rnbinom(random_n, size=random_size, prob=random_prob)
  )
}
```

Comprovació (el valor hauria de ser igual al del full de l'entregable): <ch>$\sum_{i=1}^{29} X^2 = `r format(sum(dd$v^2), scientific = FALSE)`$</ch>

## Apartat a)
**Estimeu pel mètode dels moments els paràmetres $p$ i $r$, utilitzant els dos primers moments. Per veure si l'estimació ajusta bé, dibuixeu i compareu la funció de distribució estimada $f_{mm}(x)$ amb l'empírica $F_{fr}(x)$ i contesteu les preguntes següents:**

Els moments teòrics són:

$$\begin{cases}
  \mu_1 = \mathbb{E}[X] = \frac{r(1 - p)}{p} \\
  \mu_2 = \mathbb{E}[X^2] = \mathbb{E}[X]^2 + Var(X) = \left( \frac{r(1 - p)}{p} \right)^2 + \frac{r(1 - p)}{p^2}
\end{cases}$$

Pel mètode dels moments, igualem els moments mostrals i teòrics:

$$\begin{cases}
  \hat{\mu}_1 = \frac{r(1 - p)}{p} \\
  \hat{\mu}_2 = \left( \frac{r(1 - p)}{p} \right)^2 + \frac{r(1 - p)}{p^2} = \hat{\mu}_1^2 + \frac{\hat{\mu}_1}{p}
\end{cases}$$

De la segona equació podem aïllar $p$ i trobem:

$$\hat{p}_{mm} = \frac{\hat{\mu}_1^2}{\hat{\mu}_2 - \hat{\mu}_1^2}$$

Després substituim aquest resultat a la primera equació, i aïllant $r$ trobem que:

$$\hat{r}_{mm} = \frac{\hat{\mu}_1^2}{\hat{\mu}_2 - \hat{\mu}_1^2 - \hat{\mu}_1}$$

### Pregunta 1
El valor estimat de $p$ és:

```{r}
mu1 = mean(dd$v)
#mu2 = mean(dd$v)^2 + var(dd$v) # WRONG
mu2 = sum(dd$v^2)/length(dd$v)
pmm = mu1/(mu2 - mu1^2)
pmm
```

### Pregunta 2
El valor estimat de $r$ és:

```{r}
rmm = mu1^2/(mu2 - mu1^2 - mu1)
rmm
```

### Pregunta 3
El valor estimat de $\mathbb{E}[X]$ és:
```{r}
mu1
```

### Pregunta 4
El valor estimat de $Var(X)$ és:
```{r}
# Although we could also calculate varx as |mu2 - mu1^2|, it is much cleaner to
# do it this way
#varx = var(dd$v) # WRONG
#varx
mu2 - mu1^2
```

### Pregunta 5
El valor de $\max(|F_{mm}(x) - F_{fr}(x)|)$ és:

```{r}
pemp_x = function(x) {
  return(pemp(x, dd$v, discrete=TRUE))
}

pnbinom_a = function(x) {
  return(pnbinom(x, size=rmm, prob=pmm))
}

graph_max = max(dd$v) + 5

curve(pemp_x, col="red", from=0, to=graph_max, n=3333)
curve(pnbinom_a, col="blue", from=0, to=graph_max, n=3333, add=T)
legend(1, 1, legend=c("pemp(x)", "pnbinom_mm(x)"), col=c("red", "blue"), lty=c(1, 1))
```

```{r}
curve(abs(pemp_x(x) - pnbinom_a(x)), from=0, to=graph_max, ylab="abs(pemp(x) - pnbinom_mm(x))", n=3333)
xeq = seq(-1, graph_max, 0.5)
max(abs(pemp_x(xeq) - pnbinom_a(xeq)))
```

## Apartat b)
**Estimeu pel mètode de la màxima versemblança els paràmetres $p$ i $r$. Per veure si l'estimació s'ajusta bé, dibuixeu i compareu la funció de distribució estimada $F_mm(x)$ amb l'empírica $F_{fr}(x)$ i contesteu les preguntes següents:**

La funció de probabilitat de la distribució binomial negativa és:

$$P(X = x) = \binom{r + x - 1}{x} p^r (1 - p)^x$$

Per tant, la funció de versemblança és:

$$L(p, r) \equiv L(p, r | \underset{\sim}{X}) = \prod_{i=1}^n \frac{\Gamma(r + x_i)}{\Gamma(r) \, x_i!} p^r (1 - p)^{x_i}$$

L'estimador de màxima versemblança és el que maximitza $L(p, r)$, i ho fa si i només si maximitza la funció de log-versemblança:

$$l(p, r) := \log(L(p, r)) = \sum_{i = 1}^n \left[ \log(\Gamma(r + x_i)) - \log(\Gamma(r)) - \log(x_i!) + r \log(p) + x_i \log(1 - p) \right]$$

Si derivem la funció de la log-versemblança i igualem les derivades a 0 per trobar màxims o mínims locals, trobem les següents equacions:

$$\begin{cases}
  \displaystyle 0 = \frac{\partial l}{\partial r} = \sum_{i = 1}^n \psi(r + x_i) + n (\log(p) - \psi(r)) \\
  \displaystyle 0 = \frac{\partial l}{\partial p} = n \left( \frac{r}{p} - \frac{\bar{X}}{1 - p} \right) \\
\end{cases}$$

on $\psi(x) = \frac{\partial}{\partial x} \log(\Gamma(x))$ és la funció digamma.

Reordenant la segona equació tenim:

$$p = \frac{r}{\bar{X} + r}$$

I si substituïm aquesta darrera expressió a la primera equació, obtenim:

$$0 = \sum_{i = 1}^n \psi(r + x_i) + n \left( \log\left( \frac{r}{\bar{X} + r} \right) - \psi(r) \right) =: F(r)$$

Aquí hem seguit el procediment conegut com a "profiling": trobar la solució explícita per un dels paràmetres en funció dels altres, i substituir aquesta expressió a les altres equacions per rebaixar la dimensió.

Com ens ha quedat una equació amb una incògnita, trobarem els zeros d'aquesta darrera funció per tal de trobar el valor de $r$, i mitjançant l'anterior expressió trobarem el valor de $p$.

La gràfica de la funció anterior $F(r)$ és:

```{r}
fr_scalar = function (r, x) {
  return(sum(psigamma(x + r)) + length(x)*(log(r/(mean(x) + r)) - psigamma(r)))
}

fr_list = function (r, x) {
  val = lapply(r, fr_scalar, x)
  return(val)
}

x_l = seq(0.01, 300, length.out=3333)
y_l = fr_list(x_l, dd$v)

plot(x_l, y_l, xlim=c(0.01, 125), ylim=c(-2, 2), ylab="F(r | x)", type="l")
```

### Pregunta 6
El valor estimat de $p$ és:

```{r}
fr_uniparametric = function (r) {
  return(fr_scalar(r, dd$v))
}

r_ml_result = uniroot(fr_uniparametric, lower=0.01, upper=1e5, tol=1e-10)
r_ml = r_ml_result$root
p_ml = r_ml/(mu1 + r_ml)
p_ml
```

### Pregunta 7
El valor estimat de $r$ és:

```{r}
r_ml
```

<span style="color: #aaa;">L'error estimat és de `r r_ml_result$estim.prec`</span>

### Pregunta 8
El valor estimat de $\mathbb{E}(X)$ és:
```{r}
mu1_ml = (r_ml*(1 - p_ml))/p_ml
mu1_ml
```

### Pregunta 9
El valor estimat de $Var(X)$ és:

```{r}
var_ml = (r_ml*(1 - p_ml))/(p_ml^2)
var_ml
```

### Pregunta 10
El valor de $\max(|F_{mv}(x) - F_{fr}(x)|)$ és:

```{r}
# |pemp_x| definida a la pregunta 5

pnbinom_b = function(x) {
  return(pnbinom(x, size=r_ml, prob=p_ml))
}

curve(pemp_x, col="red", from=0, to=graph_max, n=3333)
curve(pnbinom_b, col="blue", from=0, to=graph_max, n=3333, add=T)
legend(1, 1, legend=c("pemp(x)", "pnbinom_ml(x)"), col=c("red", "blue"), lty=c(1, 1))
```

```{r}
curve(abs(pemp_x(x) - pnbinom_b(x)), from=0, to=graph_max, ylab="abs(pemp(x) - pnbinom_ml(x))", n=3333)
# |xeq| definida a la pregunta 5
max(abs(pemp_x(xeq) - pnbinom_b(xeq)))
```

### Pregunta 11
Assumint que els valors dels paràmetres de $X$ fossin els valors estimats per màxima versemblança, trobeu $[a, b]$, la regió de predicció de $X$ (amb probabilitat ≥ <ch>`r prob_11`</ch> i dues cues equiprobables):

**Resolució:**

Volum una regió central amb probabilitat més gran o igual a `r prob_11` i per tant dues cues amb probabilitat menor o igual a $\alpha := \frac{1 - `r prob_11`}{2}$. A més a més, volem que les cues siguin el més gran possibles dins del rang de probabilitats que poden prendre.

```{r}
alpha_11 = (1-prob_11)/2
#a_11 = qnbinom(alpha_11, size=r_ml, prob=p_ml) + 1

a_prov = -1
for (i in 0:10000) {
  if (pnbinom(i, size=r_ml, prob=p_ml) <= alpha_11) {
    a_prov = i + 1
  } else {
    break
  }
}
a_11 = a_prov
b_11 = qnbinom(1 - alpha_11, size=r_ml, prob=p_ml)
```

> $[`r a_11`, `r b_11`]$

### Pregunta 12
De la regió de l'apartat 11) calculeu $Pr(X \in [a, b])$:

```{r}
p_12 = pnbinom(b_11, size=r_ml, prob=p_ml) - pnbinom(a_11 - 1, size=r_ml, prob=p_ml)
p_12
```


## Apartat c)
**En el supòsit que considerem que $p$ és una v.a. amb funció de densitat proporcional a $f(p) \propto p^a (1 - p)^b$, amb <ch>$a = `r a_c`$</ch> i <ch>$b = `r b_c`$</ch> amb el valor de $r$ conegut i igual a <ch>$`r r_c`$</ch>. Un cop coneguda la vostra mostra, calculeu la funció de densitat a posteriori de $p$, dibuixeu-la i contesteu les preguntes següents:**

### Pregunta 13
L'esperança de la distribució a posteriori de $p$ és:

```{r}
# Paràmetres de la distribució a posteriori:
alpha_post = a_c + length(dd$v)*r_c + 1
beta_post = b_c + length(dd$v)*mean(dd$v) + 1

esp_13 = (alpha_post)/(alpha_post + beta_post)
esp_13
```

### Pregunta 14
La variància de la distribució a posteriori de $p$ és:

```{r}
var_14 = (alpha_post*beta_post)/((alpha_post + beta_post)^2*(alpha_post + beta_post + 1))
var_14
```

### Pregunta 15
A posteriori, l'esperança de $X$ és:

```{r}
esp_15 = r_c*beta_post/(alpha_post - 1)
esp_15
```

### Pregunta 16
A posteriori, la variància de $X$ és:

```{r}
#var_16 = r_c*(r_c + 1)*(alpha_post + beta_post - 1)*beta_post/((alpha_post - 1)*(alpha_post - 2)) - esp_15^2
#var_16

r_c*(alpha_post + beta_post - 1)*beta_post/((alpha_post - 1)*(alpha_post - 2)) + r_c^2*beta_post*(beta_post + 1)/((alpha_post - 1)*(alpha_post - 2)) - r_c^2*beta_post^2/((alpha_post - 1)^2)
```

<!-- Custom styles -->
<style>
/* ch element (|ch| stands for check) */
ch {
  color: red;
}
</style>