A continuació es presenten els resultats d’aquesta pràctica amb els càlculs en R integrats en les explicacions de la solució de forma clara i llegible. En tots els exercicis tenim el conjunt de Dades 1 i venen importats cadascun. Les dades es diuen: dades1, taula_1, Dat_1 en ordre de l’exercici i segons els noms donats dels fitxers de dades.
Sigui \(X_1, \ldots, X_n\) una m.a.s. d’una v.a. \(X\) que segueix la següent distribució: \[ f(x; θ) = \frac{1}{\theta x^{\frac{1}{\theta}+1}} \qquad x>1, \quad \theta > 0 \\ E(X) = \frac{1}{1 - \theta} \qquad E(ln(X)) = \theta \quad Var(ln(X)) = \theta^2 \]
Suposarem que \(0 < θ < 1\). Donada una m.a.s \(X_1, \ldots , X_n\) de grandària \(n\) de \(X\)
Partint de l’esperança de \(X\), aïllem \(\theta\):
\[ E(X) = \frac{1}{1 - \theta} \rightarrow 1 - \theta = \frac{1}{E(X)} \xrightarrow{\text{plug-in}} \hat{\theta}_1 = 1 - \frac{1}{\overline{x}} \]
Establim la màxima versemblança amb la funció:
\[ L(x;\theta) = \prod_{i=1}^{n} f(x_i;\theta) = \prod_{i=1}^{n} \frac{1}{\theta x_i^{\frac{1}{\theta}+1}} = {\theta^{-n} \prod_{i=1}^{n}{x_i}^{-({\frac{1}{\theta} + 1})}} \]
Apliquem el logaritme neperià (i les seves propietats) per acabar de desenvolupar la funció likelihood (utilitzem que \(\ln L = l\)).
\[ \ln {L(x;\theta)} = l(x;\theta) = -n\ln{\theta} - (\frac{1}{\theta}+1)\sum_{i = 1}^{n}\ln{x_i} \] Derivant respecte \(\theta\) i igual·lant a 0, trobem el màxim versemblant \(\hat{\theta}_2\):
\[ \frac{\partial}{\partial \theta} l(x; \theta) = -\frac{n}{\theta} + \frac{1}{\theta^2}\sum_{i = 1}^{n}\ln{x_i} = 0 \] \[ -n\theta + \sum_{i = 1}^{n}\ln{x_i} = 0 \xrightarrow{\text{plug-in}} \hat{\theta}_2 = \frac{1}{n}\sum_{i = 1}^{n}\ln{x_i} \]
Calculem la fórmula de la matriu d’informació de Fischer. \[ I(\theta) = -E[\frac{\partial^2}{\partial\theta^2}\ln L(x;\theta)] \] Al tenir la funció \(L(x;\theta)\) un sol paràmetre, sabem que això donarà una matriu 1x1 d’un únic valor. Retrobem la primera derivada de l’anterior apartat (sense igual·lar a 0 aquest cop): \[ \frac{\partial}{\partial \theta} l(x; \theta) = -\frac{n}{\theta} + \frac{1}{\theta^2}\sum_{i = 1}^{n}\ln{x_i} \longrightarrow \frac{\partial^2}{\partial \theta^2} l(x; \theta) = \frac{n}{\theta^2}-\frac{2}{\theta^3}\sum_{i = 1}^{n}\ln x_i \] Acabem amb el resultat fent l’esperança (en negatiu) d’aquesta segona derivada: \[ I(\theta) = -E[\frac{n}{\theta^2}]+E[\frac{2}{\theta^3}\sum_{i = 1}^{n}\ln x_i] = -\frac{n}{\theta^2} + \frac{2}{\theta^3}\sum_{i=1}^{n}E[\ln {x_i}] \] L’esperança de les constants \(n\) i \(\theta\) són les pròpies constants. Per l’enunciat sabem també que \(E(\ln{x} = \theta)\), per tant podem acabar de resoldre. \[ I(\theta) = -\frac{n}{\theta^2} - +\frac{2}{\theta^3}\sum_{i=1}^{n}\theta = -\frac{n}{\theta^2} + \frac{2n}{\theta^2} = \frac{n}{\theta^2} \]
Quan la mostra \(n\) és suficientment gran, sabem que qualsevol paràmetre de màxima versemblança \(\hat{\theta}_2\) tendeix a la següent normal:
\[ \hat{\theta}_2 \sim N(\theta, I(\theta)^{-1}) = N(\theta, \frac{\theta^2}{n}) \]
Importem les dades del fitxer donat Dades1.csv (com a Dades1) i les integrem en aquest projecte, per així no haver d’adjuntar el fitxer csv també.
Estimació puntual de \(\theta\): Com l’estimador de màxima versemblança no té biaix ja que: \[ E(\hat{\theta}_2) = \frac{1}{n} \sum_{i = 1}^{n}E(\ln x_i) = \frac{1}{n}\sum_{i = 1}^{n}\theta = \theta \] Estimem puntualment fent \(\theta = \hat{\theta}_2\) i calculant el resultat amb R:
n <- nrow(Dades1)
theta_hat_2 <- sum(log(Dades1$x)) / n
paste("L'estimador de θ és", round(theta_hat_2, 3), "en la mostra donada de", n, "valors.")## [1] "L'estimador de θ és 0.32 en la mostra donada de 100 valors."Per \(IC_{95\%}\) de \(\theta\): Utilitzant la distribució de la màxima versemblança de \(\theta\) \[ \hat{\theta}_2 \sim N(\theta, \frac{\theta^2}{n}) \] Sabem, per la fórmula dels intervals de confiança que: \[ IC_{95\%} (\theta) = \hat{\theta}_2 \pm SE =\hat{\theta}_2 \pm z_{0.975}\frac{\hat{\theta}_2}{\sqrt{n}} \] La calculem amb R:
SE <- theta_hat_2/sqrt(n)
IC_95 <- c(theta_hat_2 - qnorm(0.975)*SE, theta_hat_2 + qnorm(0.975)*SE)
paste("L'interval de confiança del 95% de θ és", "[", round(IC_95[1],3), ",", round(IC_95[2],3), "]")## [1] "L'interval de confiança del 95% de θ és [ 0.257 , 0.383 ]"En un estudi de resultats de salts d’esquí de 61 atletes a Lilehammer es volen comparar les distàncies mitjanes (en metres) que salten homes (sexe=1) i dones (sexe=0). Fent l’estadística descriptiva, s’obtenen els resultats que apareixen a la taula 1:
## Sexe N Mean St.Dev Mínim Mediana Màxim
## 1 0 30 93.93 3.69 89.0 93.5 103
## 2 1 31 96.64 3.39 91.5 96.5 103Volem comprovar si les distàncies mitjanes de salt d’homes i dones són iguals o diferents, això planteja les següents hipòtesis:
\[ \begin{cases} H_0: \mu_1 = \mu_0 \\ H_1: \mu_1 \neq \mu_0 \end{cases} \]
Estudiem la diferència entre les distribucions de les mitjanes de salts (suposant homecedasticitat):
\[ \overline{X}_0 \sim N(\mu_o, \frac{\sigma^2}{n_0}), \quad \overline{X}_1 \sim N(\mu_1, \frac{\sigma^2}{n_1}) \]
\[ \overline{X}_1 - \overline{X}_0 \sim N(\mu_1 - \mu_0, \frac{\sigma^2}{n_0} + \frac{\sigma^2}{n_1}) \]
Al ser les variàncies iguals, però desconegudes estimem la variança a partir de l’estimació pooled calculada amb la següent fórmula.
\[ \sigma = S_p = \sqrt{\frac{(n_1 - 1)s_1^2 + (n_0 - 1)s_0^2}{n_1 + n_0 - 2}} \]
Tipificant, i sota \(H_0: \mu_1 = \mu_0\) l’estadístic de prova ens queda:
\[ t = \frac{\bar{X}_1 - \bar{X}_0}{S_p \sqrt{\frac{1}{n_1} + \frac{1}{n_0}}} \]
on \(S_p\) és la desviació típica combinada (pooled standard deviation):
Sota \(H_0: \mu_1 = \mu_0\), l’estadístic \(t\) segueix una distribució t de Student amb:
\[ gl = n_1 + n_0 - 2 \]
graus de llibertat.
Per saber-ho fem el test sota la t-Student, i calculem els p-valors per dictaminar si rebutgem \(H_0\) o no.
n0 <- taula1$N[taula1$Sexe == 0]; n1 <- taula1$N[taula1$Sexe == 1]
mean0 <- taula1$Mean[taula1$Sexe == 0]; mean1 <- taula1$Mean[taula1$Sexe == 1]
sd0 <- taula1$St.Dev[taula1$Sexe == 0]; sd1 <- taula1$St.Dev[taula1$Sexe == 1]
# Desviació típica combinada 'pooled'
sp <- sqrt(((n1 - 1)*sd1^2 + (n0 - 1)*sd0^2) / (n1 + n0 - 2))
# Estadístic t
t_stat <- (mean1 - mean0) / (sp * sqrt(1/n1 + 1/n0))
df <- n1 + n0 - 2
# p-valor (bilateral)
p_valor <- 2*(1 - pt(t_stat, df))
if (p_valor < 0.05) {
paste("Rebutgem H0 per p-valor de", round(p_valor, 5), "menor a α = 0.05")
} else {
paste("No podem rebutjar H0 per p-valor de", round(p_valor,5), "major o igual a α = 0.05")}## [1] "Rebutgem H0 per p-valor de 0.00408 menor a α = 0.05"Resultats obtinguts:
A nivell de significació \(\alpha = 0.05\), rebutgem \(H_0\) ja que \(p < 0.05\).
Com que rebutgem \(H_0\), hi ha evidència suficient per afirmar que \(H_1: \mu_1 \neq \mu_0\) és més plausible. En altres paraules, hi ha evidència que els homes salten, de mitjana, una distància diferent que les dones.
Per comprovar si és raonable assumir igualtat de variàncies, fem el test de Fisher, assumint un valor de \(\alpha = 0.05\) plasmant les següents hipòtesis:
\[ \begin{cases} H_0: \frac{\sigma_1}{\sigma_0} = 1 \\ H_1: \frac{\sigma_1}{\sigma_0} \neq 1 \end{cases} \]
Per definició de Fischer sabem que: \[ \frac{\frac{(n_0 - 1)}{(n_0 - 1)} \frac{s_0}{\sigma_0}}{\frac{(n_1 - 1)}{(n_1-1)}\frac{s_1}{\sigma_1}} = \frac{s_0 \sigma_1}{s_1\sigma_0} \sim F_{n_0-1, n_1 - 1} \] Sota \(H_0 = \frac{\sigma_1}{\sigma_0} = 1\), trobem l’estadístic de prova:
\[ F_0 = \frac{s_1^2}{s_0^2} \]
i rebutgem \(H_0\) quan l’estadístic de prova cau a les cues de la \(F_{n_0-1, n_1 - 1}\). Calculem el p-valor amb la següent fórmula:
\[ p\text{-valor} = 2 \cdot \min \Big( P(F_{n_0 -1, n_1 - 1} \le F_{0}), \; P(F_{n_0 -1, n_1 - 1} \ge F_0) \Big) \]
F_stat <- (sd1^2) / (sd0^2)
p_valor_F <- 2 * min(pf(F_stat, n1 - 1, n0 - 1), 1 - pf(F_stat, n1 - 1, n0 - 1))
if (p_valor_F < 0.05) {
paste("Rebutgem H0 per p-valor de", round(p_valor_F,3), "menor a α = 0.05")
} else {
paste("No podem rebutjar H0 per p-valor de", round(p_valor_F,3), "major o igual a α = 0.05")}## [1] "No podem rebutjar H0 per p-valor de 0.646 major o igual a α = 0.05"Resultat:
El \(p\)-valor és 0.646.
Com que és molt superior a 0.05, no es rebutja
la hipòtesi d’igualtat de variàncies, i l’assumpció
d’homocedasticitat és raonable.
Plantegem les dues hipòtesis possibles \[ \begin{cases} H_0: \frac{\mu_1} = {\mu_0} \\ H_1: \frac{\mu_1} > {\mu_0} \end{cases} \]
Ens adonem que aquest test és igual al test bilateral fet a l’apartat a), però ara hem de tractar l’unilateral cap a la dreta. Per tant reprenem l’estadístic de prova: \(T = 2.989\) i calculem el seu \(p\)-valor.
paste("El p-valor és de:", round(1 - pt(2.989, 59), 5))## [1] "El p-valor és de: 0.00204"Interpretació: Com el \(p\)-valor és menor a 0.05, cau a la zona de rebuig de \(H_0\) i podem afirmar que els homes en mitjana salten més que les dones amb un 95% de nivell de confiança.
Per afirmar si el diàmetre de les peces és de 5 cm, plantegem les següents hipòtesis:
\[ \begin{cases} H_0: \mu_0 = 5 \\ H_1: \mu_0 \neq 5 \end{cases} \]
Assumint normalitat, sabem que l’estadístic de prova és la t-Student amb n - 1 graus de llibertat, ja que estimem la variança \(\sigma = s\).
\[ T = \frac{\bar X - 5}{S / \sqrt{n}} \sim t_{n-1} \text{ sota } H_0 \]
Fem el test i obtenim tan el \(p\)-valor com l’IC del 95%.
t_res <- t.test(Dat1, mu = 5, alternative = "two.sided")
t_res##
## One Sample t-test
##
## data: Dat1
## t = 2.163, df = 49, p-value = 0.03545
## alternative hypothesis: true mean is not equal to 5
## 95 percent confidence interval:
## 5.000683 5.018585
## sample estimates:
## mean of x
## 5.009634Interpretació: com el p-valor < 0.05, rebutgem (\(H_0\)) i concluïm que la mitjana difereix de 5 cm. A més, podem apreciar com el valor 5 cm no es trobat en el 95%, i per tant, la mitjana no és de 5 cm.
Estimador puntual de la variància: \(\hat\sigma^2 = s^2\) Interval de confiança del 95% per \(\sigma^2\):
\[ IC_{95\%}(\sigma^2) = \left( \frac{(n-1)s^2}{\chi^2_{1-\alpha/2,\,n-1}},\; \frac{(n-1)s^2}{\chi^2_{\alpha/2,\,n-1}} \right) \]
alpha <- 0.05
var_x <- var(Dat1)
n <- length(Dat1)
IC_sigma <- c((n-1)*var_x / qchisq(1-alpha/2, n-1),
(n-1)*var_x / qchisq(alpha/2, n-1))
paste("L'interval del 95% de confiança de la variança és [", round(IC_sigma[1], 6), round(IC_sigma[2], 6), "]")## [1] "L'interval del 95% de confiança de la variança és [ 0.000692 0.00154 ]"Volem trobar la probabilitat de que el diàmetre de les pexes caigui fora de l’interval ([5 - 2·0.03, 5 + 2·0.03] = [4.94, 5.06]). Suposant que \(X\sim N(\mu=\bar X,\sigma= s)\) (utilitzem les estimacions mostral):
fita_inferior <- 4.94; fita_superior <- 5.06
mitja_x <- mean(Dat1)
sd_x <- sqrt(var_x)
p_defectuosa <- pnorm(fita_inferior, mitja_x, sd_x) + (1 - pnorm(fita_superior, mitja_x, sd_x))
paste("La probabilitat de ser defectuosa és de",round(p_defectuosa,4))## [1] "La probabilitat de ser defectuosa és de 0.0684"Com la probabilitat de ser defectuosa, que hem acabat de calcular, és molt major a 0.03, intuïm que no podem considerar que sigui inferior al 3%. No obstant, com això és una pràctica d’Estadística hem de verificar-ho estadísticament! Formulem el test (unidireccional) on volem veure si és la probabilitat de defectuosa és inferior a 3% (hipòtesi alternativa):
\[ \begin{cases} H_0: p_{0} = 0.03 \\ H_1: p_{0} < 0.03 \end{cases} \]
Utilitzem el valor estimat de (p = p_{defectuosa}), calculat a c) i fem una aproximació de la distribució binomial de les peces defectuoses a una normal de paràmetres:
\[ \hat{p} \sim N(p_{0}, \frac{p_0(1-p_0)}{n}) \] Finalment, calculem els p-valors, fent el test sota aquesta distribució:
p0 <- 0.03
p_hat <- p_defectuosa
z_stat <- (p_hat - p0) / sqrt(p0*(1-p0)/n)
p_valor_test <- pnorm(z_stat) # unidireccional esquerra
paste("Test amb p-valor:", p_valor_test)## [1] "Test amb p-valor: 0.944319142076625"Interpretació: donat que el p-valor > 0.05, no podem afirmar que la probabilitat de trobar una peça defectuosa sigui inferior al 3%.
Establim dues variables aleatòries que mesuren les mitjanes dels diàmetres de les peces, assumint que la variança és la mateixa (ja que mesurem la mateixa distribució de diàmetre de peces), i imaginem que una es descentra 0.5 cm’s de l’altre
\[ \overline{D_1} \sim N(\mu_1, \frac{\sigma^2}{n}), \qquad \overline{D_2} \sim N(\mu_1 + 0.5, \frac{\sigma^2}{n}) \\ \overline{D_1}-\overline{D_2} \sim N(0.5, 2\frac{\sigma^2}{n}) \]
Sigui quina sigui la \(n\), assegurem que l’error de tipus I sigui \(\alpha = 0.4\). Plantegem que a \(\overline{D_1}\) rebutgem \(H_0\) quan sobrepassa el valor crític:
\[ z_{1-\frac{\alpha}{2}} = z_{0.8} = \frac{\overline{D_1}-\overline{D_2}-\mu_1}{\sigma\sqrt{\frac{2}{n}}} \] I no rebutgem quan (fixant el valor de \(\beta = 0.4\):
\[ z_{\beta} = z_{0.4} = \frac{\overline{D_1}-\overline{D_2}-(\mu_1+k)}{\sigma\sqrt{\frac{2}{n}}} \]
Restem ambdues equacions per trobar el valor de \(n\), i estimant la variança per la variança mostral i la mitjana per la mitja mostral.
\[ z_{0.8} - z_{0.4} = \frac{0.5}{s\sqrt{\frac{2}{n}}} \\ n = (\frac{(z_{0.8}-z_{0.4})s\sqrt{2}}{0.5})^2 \]
n <- (((qnorm(0.8)-qnorm(0.4))*sd_x*sqrt(2))/0.5)^2
paste("Ha de tenir com a mínim grandària", round(n, 4))## [1] "Ha de tenir com a mínim grandària 0.0095"Resultat: Necessitariem n = 0.0095 \(\approx\) 1 mostra per detectar diferències de 0.5 cm, amb error de tipus I i II de 0.4.
L’enunciat proposa un model binomial de 50 peces el qual aproximem com una normal de paràmetres:
\[ X \sim N(n p_{def}, n{p_{def}(1-p_{def})}) \] Per tant, per trobar la probabilitat de tenir 3 o més peces defectuoses, és la mateixa de no tenir 2 o menys peces defectuoses, que calculat amb R dona:
p_3_o_mes <- 1- pnorm(2, mean = 50*p_defectuosa, sd = sqrt(50*p_defectuosa*(1-p_defectuosa)))
paste("La probabilitat de que hi hagi 3 defectuoses o més en una caixa de 50 peces és de ", round(p_3_o_mes, 3))## [1] "La probabilitat de que hi hagi 3 defectuoses o més en una caixa de 50 peces és de 0.787"