\documentclass[12pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[spanish]{babel}
\usepackage[margin=2cm]{geometry}
\begin{document}
\title{\bf Segundo parcial}
\date{1 de Junio de 2017}
\maketitle
\begin{enumerate}
\item (2 Puntos) Explica detalladamente el test de Mann-Whitney:
  cuáles son las condiciones de aplicación, el estadístico de
  contraste y la región crítica. Justifica cómo se calcula la
  distribución del estadístico bajo la hipótesis nula si $n_1=n_2=5$.

\item (2 Puntos) En un análisis ANOVA con $5$ niveles y $10$
  observaciones por nivel explica cómo se define el estadístico SCF.
  ¿Qué estima y cuál es su distribución bajo la hipótesis nula?
  ¿Por qué no se puede utilizar sólo este estadístico para el contraste
  ANOVA?

\item (3 Puntos) La variación en milímetros de la anchura, $X$, de
  cierta pieza respecto al patrón de fabricación puede considerarse
  una variable aleatoria con densidad
  $$f(x)= \frac{|x|}{\theta^2} I( |x|\le\theta) $$
  A partir de la información muestral siguiente
  \begin{flushleft}
  {\tt (-1.85, -1.37, 2.28,  2.09, -2.23,  1.91, -0.34, -0.63, -0.37,  1.71, 2.25,  2.16, 1.90, -2.25,
  2.10, 2.47, -2.27,  1.95, -1.14, -1.70)}
\end{flushleft}
  responde a las siguientes preguntas:
\begin{enumerate}

\item Contrasta las hipótesis $H_0:\theta \le 2.5$ frente
  $H_1:\theta > 2.5$ al nivel de significación $0.05$ basándote en el
  estadístico de la razón de verosimilitudes.

\item Contrasta si los datos anteriores provienen de la familia de
  distribuciones especificada en el enunciado, aplicando un test
  chi-cuadrado.

\end{enumerate}

\item (3 Puntos) Para comparar el nivel de contaminación en cuatro
  plantas de energía diferentes (Plantas: 1,2,3,4) se realizó un
  experimento en el que se medía el nivel máximo de SO$_2$ diario y el
  grado de contaminación ambiental (Bajo, Medio, Alto, Muy Alto). Los
  resultados obtenidos están en el fichero: \textit{contamina.rda}

  \begin{enumerate}
  \item Analiza si existen diferencias de SO$_2$ según el tipo de
    planta.
  \item Contrasta si el grado de contaminación es similar en todas las
    plantas.
  \item Estudia si existe relación entre el nivel máximo de SO$_2$ y
    el grado de contaminación ambiental.
  \item Comenta las conclusiones que obtienes del experimento.
\end{enumerate}

\end{enumerate}
JUSTIFICA Y COMENTA TODOS LOS RESULTADOS.

\newpage

\section*{Solución}

\subsection*{Ejercicio 3}
\subsubsection*{Apartado a)}
Se quieren contrastar las hipótesis $ H_0:\theta \le 2.5$ frente a
$H_1:\theta > 2.5$ al nivel de significación $\alpha=0.05$. Para
buscar una región crítica se pueden utilizar distintas alternativas,
como por ejemplo:
\begin{itemize}
\item Calcular un buen estimador del parámetro y decidir qué evidencia
  es necesaria para rechazar $H_0$.
\item Utilizar el test de la razón de verosimilitudes.
\item Por ser un test unilateral se puede intentar emplear el criterio
  de Karlin-Rubin.
\end{itemize}

\newcommand{\maxabs}{{|X|_{(n)}}}
En primer lugar se calcula el estimador máximo-verosímil de $\theta$:
$$
L(x_1, ..., x_n, \theta) = 
\prod_{i=1}^n  \frac{|x_i|}{\theta^2} \, I(|x_i|\le \theta ) =
\prod_{i=1}^n  \frac{|x_i|}{\theta^2} \, I\left(\max_i |x_i|\le \theta \right) =$$
$$= \left(\frac{1}{\theta}\right)^{2n} \, \prod_{i=1}^n  |x_i| \, I(\maxabs\le \theta),  \quad 
\mbox { con } \maxabs =\max_i |x_i|
$$

La función de verosimiltud alcanza su máximo en el valor más pequeño
del parámetro que cumpla la condición $\maxabs\le\theta$. Por tanto el
estimador máximo-verosímil es:
$$ \mbox{EMV}(\theta)=\hat{\theta} = \maxabs $$

Teniendo en cuenta cómo son las hipótesis, lo que
debería ocurrir bajo $ H_0:\theta \le 2.5$ es que el EMV tome valores
más pequeños que bajo $H_1: \theta > 2.5$. Por tanto una región
crítica razonable es $$\mbox{RC}= \{ \maxabs > k \}$$

La distribución de $X$ verifica la propiedad de razón de verosimilitud
monótona respecto a $\maxabs$. Sean $\theta_2>\theta_1$; entonces
$$
\frac{L(\vec{x},\theta_2)}{L(\vec{x},\theta_1)}=  
\left(\frac{\theta_1}{\theta_2}\right)^n
\frac{I(\maxabs \le \theta_2)}{I(\maxabs \le \theta_1)}
 = \left\{
\begin{array}{cl}
  \left(\frac{\theta_1}{\theta_2}\right)^n &  \maxabs \le \theta_1\\
  \infty & \theta_1 < \maxabs \le \theta_2 \\
\end{array}  \right.
$$

Por tanto se verifican las condiciones para aplicar el teorema de
Karlin-Rubin y la región crítica del test es de la forma:
$$\mbox{RC} = \{\maxabs > h\}$$

Empleando el método de la razón de verosimilitudes se obtiene la
misma región crítica.  Téngase en cuenta que si se expresa la RC en
función de $\Lambda$ surge la necesidad de determinar un procedimiento
aleatorizado, dado que la razón de verosimilitudes toma valor
constante en el intervalo de $\maxabs$ que produce verosimilitudes
no nulas.

La siguiente etapa es calcular la distribución del estadístico de
contraste, bajo $H_0$, y determinar el valor de $k$ para
$\alpha=0.05$. Para ello hay que calcular la distribución de
$Y = |X|$.

$$
P(Y \le y) = P(|X| \le y) = P( -y \le X \le y) 
= \int_{-y}^{y} \frac{|x|}{\theta^2}\, dx
= \int_{0}^{y} 2\, \frac{x}{\theta^2}\, dx
= \frac{y^2}{\theta^2}
$$ 

$$ P(RC \mid H_0) = P(\maxabs > h \mid H_0)=  
1 - P(\maxabs \le h \mid H_0)=
1 - P(|X| \le h \mid H_0)^n=
$$
$$= 1- \left(\frac{h^2}{\theta^2} \right)^n 
\le 1- \left(\frac{h}{\theta_0} \right)^{2n} \le \alpha
, \quad \mbox{ con } \theta_0 = 2.5$$

Por tanto se tiene que
$\left(\frac{h}{\theta_0} \right)^{2n}= 1-\alpha \Longleftrightarrow
h= \theta_0 \, \left( 1-\alpha \right)^{1/2n}=2.5\times0.95^{1/40} = 2.4968$

Cálculos en R:
\begin{verbatim}
## Ejercicio 3
## Apartado a)
## H0: theta <= theta0; H1: theta > theta0; theta0 = 2.5 
library(RcmdrMisc)               # para numSummary
theta0 <- 2.5
alfa   <- 0.05
x      <- c(-1.85, -1.37, 2.28,  2.09, -2.23,  1.91, -0.34, -0.63, 
            -0.37,  1.71, 2.25,  2.16, 1.90, -2.25,
            2.10, 2.47, -2.27,  1.95, -1.14, -1.70)
numSummary(x)                    # descripción de los datos
n <- length(x)                   # tamaño de la muestra
theta_e <- max(abs(x))           # estimador máximo-verosímil
 
## RC = {máx |Xi| > h}           # región crítica
h <- theta0 * (1-alfa)^(1/(2*n)) # frontera de la RC
\end{verbatim}
\subsubsection*{Apartado b)}
Para usar el test chi-cuadrado de bondad de ajuste con una variable
continua, hay que discretizar la variable en categorías cuya
frecuencia esperada sea mayor o igual a cinco, bajo la hipótesis nula.

Para estimar los umbrales de las categorías se calcula la función de
distribución de $X$ sustituyendo el valor de $\theta$ por el EMV. Como
la distribución es simétrica con mediana igual a cero, sólo es
necesario hallar los puntos de corte para valores negativos.
$$ F_X(x) = -\frac{-1}{2} \, \left(\frac{x^2}{\theta^2} -1\right)=
\frac12 \left(\frac{x^2}{\theta^2} -1\right), \quad x \le 0$$

Si elegimos discretizar la variable en cuatro (=$n/5$) categorías
equiprobables, tendríamos que calcular el cuantil $x_p$ de la variable
asociado a $p=0.25$. Los puntos de corte serían entonces $-$EMV,
$x_p$, $0$, $-x_p$, EMV.

$$ p=\frac12 \left(\frac{x^2}{\theta^2} -1\right)
\Longrightarrow x_p = -\theta \sqrt{1-2*p}$$

Con la muestra del enunciado:
\begin{verbatim}
## Apartado b)
numSummary(x)                  # descripción de los datos 
hist(x)
boxplot(x)
## puntos de corte
P1     <- 0.25
pc1    <- -theta_e * sqrt(1-2*P1) # punto de corte para F(x)=0,25
cortes <- c(-2.5, pc1, 0, -pc1, 2.5)
xc     <- cut(x, cortes)
T      <- table(xc)                 # frecuencias observadas
P0     <- rep(.25, 4)               # probabilidades bajo H0
E      <- n * P0                    # frecuencias esperadas bajo H0
esta   <- sum((T-E)^2 / E)          # estadísico de Pearson
chisq.test(table(xc), p=rep(.25,4)) # usa 3 grados de libertad
## X-squared = 6.8, df = 3, p-value = 0.07855
1 - pchisq(esta, 2)                 # 2 g.l. por estimar 1 parámetro
## p-valor = 0,03337
## pero deberíamos estimar el parámetro con datos agrupados
F <- function (x, theta)            # función de distribución
  ifelse (x < 0,                    #   entre -theta y +theta
          -x^2 / (2*theta^2) + 1/2, # si -theta < x < 0
          1/2 + x^2 / (2*theta^2))  # si 0 < x < theta
logL <- function (theta)            # logverosimilitud
  sum(T * log(diff(F(c(-theta, pc1, 0, -pc1, theta), theta))))
theta_a <- optimize (logL, c(-pc1,10), maximum=TRUE) $ maximum
## estimación = 2,952
E       <- 20 * diff(F(c(-theta_a, pc1, 0, -pc1, theta_a), theta_a))
## E = (6,5 3,5 3,5 6,5) ¡ojo!
1 - pchisq (sum ((T-E)^2 / E), 2)
## p-valor = 0,0641
\end{verbatim}

\subsection*{Ejercicio 4}
Consideraremos en general un nivel de significación $\alpha=0.05$.
\begin{verbatim}
load("~/Descargas/contamina.RData")

## apartado a) SO2 según las Plantas
summary(contamina)
boxplot(SO2 ~ Planta, contamina, notch=TRUE)
by(contamina$SO2, contamina$Planta, shapiro.test)
## rechazamos gausianidad incluso con Bonferroni
kruskal.test(SO2 ~ Planta, contamina)
## p-valor = 0,022
pairwise.wilcox.test(contamina$SO2, contamina$Planta)
## diferencias sólo entre Plantas 2 y 3

## apartado b) contaminación según Plantas
Ta <- table(contamina$Planta, contamina$contamina)
Ta                         # dados los pequeños tamaños
prop.table(Ta, 1)          # no se aprecian diferencias claras
chisq.test(Ta)             # aviso: alguna esperada < 5
chisq.test(Ta, simulate.p.value=TRUE) # p-valor parecido

## apartado c) SO2 según la contaminación
boxplot(SO2 ~ contamina, contamina, notch=TRUE)
by(contamina$SO2, contamina$contamina, shapiro.test)
## rechazamos gausianidad
kruskal.test(SO2 ~ contamina, contamina)
## p-valor = 0,59
\end{verbatim}
\subsubsection*{Apartado d)}
\begin{enumerate}
\item Hay evidencias leves de que las Plantas 2 y 3 tienen distintas
  distribuciones del máximo de SO$_2$.  Sin embargo, no se detectan
  diferencias entre (por ejemplo) la 2 y la 1, ni entre la 3 y la 1.
\item No hay evidencias de que el grado de contaminación difiera entre
  Plantas.
\item No hay evidencias de que el grado de contaminación esté asociado
  al máximo de SO$_2$.
\end{enumerate}
\end{document}
