Add entregables 2, 3, 4
Change-Id: I7ef88f875fdc53abcf1744a8a79a1c0a6512105b
diff --git a/entregables/ent2/ent2.Rmd b/entregables/ent2/ent2.Rmd
new file mode 100644
index 0000000..f390a54
--- /dev/null
+++ b/entregables/ent2/ent2.Rmd
@@ -0,0 +1,357 @@
+---
+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>