| --- |
| 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> |