\documentclass[12pt]{article}
\usepackage[spanish]{babel}
\usepackage[margin=1cm]{geometry}
\usepackage{palatino}
\usepackage{color}
\newenvironment{resoluciOn}
{\begin{quotation}\sffamily\color{blue}}
  {\end{quotation}}
\makeatletter
\renewcommand{\verbatim@font}{\ttfamily\small}
\makeatother
\newcommand{\betaqq}{{\rm B}(q,q)}
\newcommand{\betaprimaqq}{{\rm B}'(q,q)}
\newcommand{\sigue}{\;\hookrightarrow\;}
\begin{document}
\pagestyle{empty}
\boldmath
\noindent{\Large \bf Inferencia Estadística --- 20 de diciembre de 2024
  --- primer parcial.}

\paragraph{Instrucciones.}
\begin{itemize}
\item Primera parte: TEORÍA. Ordenador apagado. Tienes hasta una hora
  para resolverla y entregar los folios con tu respuesta.
\item Segunda parte: PROBLEMAS. En cuanto entregues la teoría, podrás
  encender el ordenador exclusivamente para usar R y acceder al campus
  virtual. No hay una ``tarea'' en el campus virtual para subir documentos.
  Has de entregar todas tus respuestas en papel: para un gráfico, esbózalo
  a mano dejando claro lo que quieras destacar; para código fuente,
  escríbelo con la mayor precisión posible.  Has de entregar los folios en
  dos bloques separados: problema 1; problema 2.
\item Justifica todas las respuestas.
\end{itemize}
\paragraph{Teoría (4 puntos).}
\begin{enumerate}
\item A una inteligencia artificial se le planteó el siguiente problema:
  \begin{quotation}
\noindent\textbf{Pregunta:}
\textit{Sea $X$ una variable aleatoria con densidad
$f(x, \lambda) = \lambda * exp(- \lambda*(x-\lambda)) \; I(x >\lambda)$. ¿Puedes calcular el estimador máximo-verosímil de $\lambda$?} % y escribir el resultado} en latex.% Puedes darme un programa en R para calcular landa para la muestra:
%$$(6.289, 1.123, 2.031, 5.822, 3.208, 2.684, 0.961, 0.698, 0.711, 3.862)$$
\end{quotation}
Analiza si la respuesta es correcta y en caso contrario comenta los
posibles fallos:
\begin{quotation}
\noindent\textbf{Respuesta:} \it
Para calcular el estimador de máxima verosimilitud (EMV) del parámetro $\lambda$ para una variable aleatoria X con densidad:
$$f(x|\lambda) = \lambda \exp(-\lambda (x - \lambda))  \quad (x > \lambda)$$, seguimos estos pasos:
\begin{enumerate}
 \item Función de Verosimilitud
%La función de verosimilitud para una muestra $\vec{x}= (x_1, x_2, \ldots, x_n)$ es:
$$L(\lambda) = \prod_{i=1}^n \lambda \exp(-\lambda (x_i - \lambda)) = \lambda^n \exp\left(-\lambda \sum_{i=1}^n (x_i - \lambda)\right)$$

\item Logaritmo de la Función de Verosimilitud
%El logaritmo de la función de verosimilitud es:
$$\log L(\lambda) = n \log \lambda - \lambda \sum_{i=1}^n (x_i - \lambda) =n \log \lambda - \lambda \left( \sum_{i=1}^n x_i - n\lambda \right)$$
% \item Maximización de la Función de Verosimilitud
%
% Para encontrar el EMV de $\lambda$, derivamos respecto a $\lambda$ y encontramos el valor que hace que la derivada sea cero:
% $$\frac{d}{d\lambda} \log L(\lambda) = \frac{n}{\lambda} - \sum_{i=1}^n (x_i - \lambda) = 0$$
%
% Resolviendo para $\lambda$:
% $$\frac{n}{\lambda} = \sum_{i=1}^n (x_i - \lambda)$$
%
% %$$\frac{n}{\lambda} = \sum_{i=1}^n x_i - n\lambda$$
% %$$\frac{n}{\lambda} = n \bar{x} - n\lambda$$
% $$\frac{1}{\lambda} = \bar{x} - \lambda$$
% $$\lambda (\bar{x} - \lambda) = 1$$
% $$\lambda = \frac{1}{\bar{x}} + \lambda$$
% $$\lambda = \sqrt{\bar{x}}$$

\item Maximización de la Función de Verosimilitud\\
Derivamos respecto a $\lambda$ y encontramos el valor que hace que la derivada sea cero:
$$\frac{d}{d\lambda} \log L(\lambda) = \frac{n}{\lambda} -
\sum_{i=1}^n x_i + 2n\lambda = 0$$
%$$\frac{n}{\lambda} = \sum_{i=1}^n x_i - 2n\lambda$$
$$\frac{n}{\lambda} + 2n\lambda = \sum_{i=1}^n x_i  \Longleftrightarrow
\lambda(1 + 2\lambda) = \frac{1}{n} \sum_{i=1}^n x_i$$
\end{enumerate}
% \begin{verbatim}
%  # Muestra de datos
% x <- c(6.289, 1.123, 2.031, 5.822, 3.208, 2.684, 0.961, 0.698, 0.711, 3.862)
%
% # Estimador de máxima verosimilitud para lambda
% lambda_hat <- mean(x)
%
% # Resultado
% lambda_hat
% \end{verbatim}

\noindent\textbf{Pregunta:}
¿Puedes resolver la ecuación $\lambda(1 + 2\lambda) = \frac{1}{n} \sum_{i=1}^n x_i$ en $\lambda$?

\noindent\textbf{Respuesta:}
Para resolver la ecuación podemos reescribirla en términos de una ecuación cuadrática estándar.
$$ 2\lambda^2 + \lambda - \frac{1}{n} \sum_{i=1}^n x_i = 0$$
%Ahora, tenemos una ecuación cuadrática en la forma estándar $a\lambda^2 + b\lambda + c = 0$, donde:

%$a=2$, $b=1$ y $c = -\frac{1}{n} \sum_{i=1}^n x_i$

% Sustituyendo $a$, $b$ y $c$ en la fórmula general, obtenemos dos posibles soluciones para $\lambda$:
Ahora, tenemos una ecuación cuadrática en la forma estándar y obtenemos dos posibles soluciones: % para $\lambda$:

$$\lambda = \frac{-1 \pm \sqrt{1 + \frac{8}{n} \sum_{i=1}^n x_i}}{4}$$
Dado que $\lambda$ representa una cantidad física y debe ser positiva, seleccionaremos la solución positiva:
$$\lambda = \frac{-1 + \sqrt{1 + \frac{8}{n} \sum_{i=1}^n x_i}}{4}$$
\end{quotation}
% \item Sea $X$ una variable aleatoria con distribución gamma
%   $\gamma(p,a)$ de la que se extrae una muestra aleatoria simple
%   $(X_1, \ldots, X_n)$ con $n$ muy grande. Calcula la distribución
%   asintótica de $\bar{X}^2$. Recuerda que $E(X) = \frac{p}{a}$ y
%   $\mathit{Var}(X) = \frac{p}{a^2}$.
\item Explica el procedimiento de \emph{bootstrap} básico para calcular intervalos de
  confianza.
\item Define partición suficiente y minimal. Demuestra que la
  partición
  $$\vec{x} \sim \vec{y}\quad \Longleftrightarrow \quad \frac{L(\vec{x},
    \theta)}{L(\vec{y}, \theta)} \mbox{ independiente de } \theta$$
  es suficiente y minimal.
\item Define la región crítica de un contraste de hipótesis y el nivel
  de significación $\alpha$.
\end{enumerate}
  \paragraph{Problema 1 (3 puntos).}
Cierta ceca fabrica monedas de forma que la probabilidad de que salga cara si tiramos
al aire la moneda $i$-ésima es $p_i$. Es decir, distintas monedas
tienen distinta probabilidad de cara.  Cada $p_i$ puede estimarse con
mucha precisión, por lo que dada una muestra
concreta de $n$ monedas sus $p_i$ se suponen conocidas ($i=1,\dots,n$).
Las $p_i$ provienen de una población beta
$\beta(q,q)$, es decir, con ambos parámetros iguales.  Se dispone de una
realización muestral de tamaño $n=10$, a saber, $(p_1,\dots,p_{10})$ =
%set.seed(1);(x<-round(rbeta(10,25,25),3));stopifnot(all(!duplicated(x)))
\texttt{\bfseries(0.449, 0.515, 0.432, 0.526, 0.433, 0.539, 0.560, 0.546,
  0.476, 0.630)}.
\begin{resoluciOn}
  En R definimos una variable que albergue la muestra:
\begin{verbatim}
> x <- c(0.449, 0.515, 0.432, 0.526, 0.433, 0.539, 0.560, 0.546, 0.476, 0.630)
> summary(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.4320  0.4557  0.5205  0.5106  0.5443  0.6300 
\end{verbatim}
\end{resoluciOn}
\begin{itemize}
\item[a)] Halla un estadístico suficiente minimal para $q$  ($0.6$ puntos).
  \begin{resoluciOn}
    Para $0<x<1$, la densidad es $$f(x)=x^{q-1}\cdot(1-x)^{q-1}/\betaqq$$
    luego la verosimilitud es
    $$\mathcal L(\vec x, q)=\prod x_i^{q-1}\cdot\prod(1-x_i)^{q-1}/\betaqq^n$$
    y el cociente de verosimilitudes es
    $$
    \frac{\mathcal L(\vec x,q)}{\mathcal L(\vec y,q)}=
    \frac{\prod x_i^{q-1}\cdot\prod(1-x_i)^{q-1}/\betaqq^n}
    {\prod y_i^{q-1}\cdot\prod(1-y_i)^{q-1}/\betaqq^n}
    =\left(\frac{\prod x_i\cdot\prod(1-x_i)}{\prod y_i\cdot\prod(1-y_i)}\right)^{q-1}
    $$
    el cual es independiente del parámetro $q$ si y sólo si
    $$
    {\prod x_i\cdot\prod(1-x_i)}={\prod y_i\cdot\prod(1-y_i)}
    $$
    luego el estadístico suficiente minimal será
    $$T(\vec x)={\prod x_i\cdot\prod(1-x_i)}=\prod[x_i(1-x_i)]$$

    Otra forma de verlo se basa en tener en cuenta que la distribución
    beta forma parte de la familia exponencial:
    $$
    \mathcal L = \frac1{\betaqq^n}\cdot\frac{1}{\prod x_i(1-x_i)}\cdot
    e^{q\cdot\ln\prod x_i(1-x_i)}
    $$
    luego $\ln\prod x_i(1-x_i)$ es suficiente minimal y, por tanto,
    $\prod x_i(1-x_i)$ también lo es por ser una trasformación biyectiva suya.
  \end{resoluciOn}
\item[b)] Halla una estimación de $q$ por el método de máxima verosimilitud
  ($0.8$ puntos).
  \begin{resoluciOn}
    Primero intentaremos encontrar el estimador máximo-verosímil:
    \begin{eqnarray*}
    \mathcal L &=& \prod x_i^{q-1}\cdot\prod(1-x_i)^{q-1}/\betaqq^n \\
      \ln\mathcal L &=& (q-1)\sum\ln x_i+(q-1)\sum\ln(1-x_i)-n\ln \betaqq \\
      0=\frac{\partial\ln\mathcal L}{\partial q}
      &=&\sum\ln x_i+\sum\ln(1-x_i)-n\frac{\betaprimaqq}{\betaqq}
    \end{eqnarray*}
    luego no se puede despejar $q$ de forma analítica, por lo que haremos
    una maximización numérica.
    
    La verosimilitud se puede codificar en R de una de estas dos maneras:
\begin{verbatim}
> n <- length(x)
> L <- function (q) prod(x)^(q-1)*prod(1-x)^(q-1)/beta(q,q)^n
> L <- function (q) prod (dbeta (x, q, q))
\end{verbatim}
    y puede maximizarse así:
\begin{verbatim}
> optimize (L, c(0,100), maximum=TRUE)
$maximum
[1] 33.08616

$objective
[1] 965217
\end{verbatim}
    luego la estimación máximo verosímil es $33.08616$.

    Una comprobación gráfica podría hacerse así:
\begin{verbatim}
> emv <- optimize (L, c(0,100), maximum=TRUE) $ maximum
> plot (Vectorize(L), 0, 100)
> abline (v=emv)
\end{verbatim}
    donde se aprecia la unimodalidad de la verosimilitud y, para comprobar
    la adecuación de la estimación a la muestra:
\begin{verbatim}
> hist (x, prob=TRUE, xlim=c(0,1))
> curve (dbeta (x, emv, emv), 0, 1, add=TRUE)
\end{verbatim}
    donde se aprecia que, al estar los datos muy cercanos $0.5$, tiene sentido
    un exponente $q$ grande en la expresión de la densidad.

    Alternativamente podría haberse usado la logverosimilitud de alguna de estas
    formas:
\begin{verbatim}
> lnL <- function (q) (sum(log(x))+sum(log(1-x)))*(q-1)-log(beta(q,q))*n
> lnL <- function (q) sum (dbeta (x, q, q, log=TRUE))
> optimize (lnL, c(0,100), maximum=TRUE) # el lugar del máximo no cambia
$maximum
[1] 33.08617

$objective
[1] 13.78011
\end{verbatim}
  \end{resoluciOn}
\item[c)] Halla una estimación de $q$ por el método de los momentos
  ($0.8$ puntos).
  \begin{resoluciOn}
    Como $E(X)=q/2q=1/2$, hay que recurrir al momento de segundo orden.

    Si empleamos un momento relativo,
    \begin{eqnarray*}
      S^2&\approx&\mathit{Var}(X) = q^2/4q^2(2q+1)=1/4(2q+1)\\
      2q+1&\approx&1/4S^2\\
      \hat q_{\rm MM}&=&1/8S^2-1/2
    \end{eqnarray*}
    Para la muestra del enunciado:
\begin{verbatim}
> 1/(8*var(x)) - 1/2 # usando cuasivarianza
[1] 30.57357
> 1/(8*(mean(x^2)-mean(x)^2)) - 1/2 # usando varianza
[1] 34.02619
\end{verbatim}

    Si usamos un momento absoluto,
    \begin{eqnarray*}
      \overline{X^2}&\approx&\mathit{Var}(X)+E(X)^2=1/4(2q+1)+1/4\\
      2q+1&\approx&1/4(\overline{X^2}-1/4)=1/(4\overline{X^2}-1)\\
      \hat q_{\rm MM}&\approx&1/(8\overline{X^2}-2)-1/2
    \end{eqnarray*}
    En este caso, la función es discontinua y tiende a $\pm\infty$
    cuando $\overline{X^2}=2$,
    por lo que el estimador puede no ser consistente.  De hecho, en este
    caso la estimación se aleja de la máximo-verosímil:
\begin{verbatim}
> 1/(8*mean(x^2)-2) - 1/2
[1] 8.221255
\end{verbatim}
  \end{resoluciOn}
\item[d)] Halla un intervalo de confianza para $q$ a nivel $95\%$
  ($0.8$ puntos).
  \begin{resoluciOn}
    No disponemos de pivote exacto.  Podríamos pensar en pivote
    asintótico (ver más adelante) pero el tamaño muestral $n=10$
    parece demasiado pequeño (aunque la distribución es simétrica,
    lo que podría ayudar a acelerar la convergencia).
    De momento, usaremos un $t$-\emph{bootstrap} paramétrico:
\begin{verbatim}
dis <- replicate (1000,
{
  xb <- rbeta (n, emv, emv)
  Lb <- function (q) prod (dbeta (xb, q, q))
  qb <- optimize (Lb, c(0,100), maximum=TRUE) $ maximum
  sigma <- sd (replicate (100, {
      xb <- rbeta (n, qb, qb)
      Lb <- function (q) prod (dbeta (xb, q, q))
      optimize (Lb, c(0,100), maximum=TRUE) $ maximum}))
  c (qb, (qb-emv)/sigma)
})
sigma <- sd (dis[1,])
emv - sigma * quantile(dis[2,], c(.975,.025), names=FALSE)
\end{verbatim}
    Al ejecutarlo varias veces se obtienen resultados como los siguientes:
\begin{verbatim}
[1] -46.47870  67.63799
[1] -45.80360  62.64771
[1] -42.15039  65.92075
\end{verbatim}
    Esos extremos inferiores negativos no cuadran mucho, porque el parámetro $q$
    ha de ser positivo.  Explorando las estimaciones obtenidas se observa que
    hay una acumulación en el extremo superior de la búsqueda numérica:
\begin{verbatim}
> stem(dis[1,])

 The decimal point is 1 digit(s) to the right of the |

 1 | 1223334444
 1 | 5555556666666666677777777777888888888888999999999
 2 | 00000000000000111111111111111111111111111222222222222222222222222222+50
 2 | 55555555555555555555555555566666666666666666666666677777777777777777+49
 3 | 00000000000000000000000000000000001111111111111111111111111111111111+81
 3 | 55555555555555555555555556666666666666666666666666677777777777777777+29
 4 | 00000000000000000000111111111111111111111122222222222223333333333333+12
 4 | 55555555555555555666666666666666677777777777788888888888888999999999
 5 | 00000000011111111222222222222333333444444
 5 | 5555555555555666677777777778888888899999999
 6 | 000000111222222223333333333333334444
 6 | 555555566666677777788888999999
 7 | 0000001111222233334444
 7 | 6677778889999
 8 | 1223444
 8 | 556667899
 9 | 11112223334
 9 | 56789
10 | 0000000000000000000000000000000000
\end{verbatim}
    Si se repite el mismo código pero ampliando a 1000 el extremo
    superior de la búsqueda numérica, ya no se obtienen extremos
    negativos.  Ampliaremos también a 10.000 el número de repeticiones
    para mejorar la precisión de la estimación \emph{bootstrap}:
\begin{verbatim}
> dis <- replicate (10000,
+ {
+   xb <- rbeta (n, emv, emv)
+   Lb <- function (q) prod (dbeta (xb, q, q))
+   qb <- optimize (Lb, c(0,1000), maximum=TRUE) $ maximum
+   sigma <- sd (replicate (100, {
+       xb <- rbeta (n, qb, qb)
+       Lb <- function (q) prod (dbeta (xb, q, q))
+       optimize (Lb, c(0,1000), maximum=TRUE) $ maximum}))
+   c (qb, (qb-emv)/sigma)
+ })
> sigma <- sd (dis[1,])
> emv - sigma * quantile(dis[2,], c(.975,.025), names=FALSE)

[1] 10.53789 69.30451
\end{verbatim}
    Ahora el comportamiento de las estimaciones parece correcto:
\begin{verbatim}
> stem(dis[1,])

  The decimal point is 1 digit(s) to the right of the |

 0 | 99000111111111222222222222222223333333333333333333333333333333444444+712
 2 | 00000000000000000000000000000000000000000000000000000000000000000000+5155
 4 | 00000000000000000000000000000000000000000000000000000000000000000000+2545
 6 | 00000000000000000000000000000000000000000000000000000000000001111111+737
 8 | 00000000000000000000000000000111111111111111111111122222222222222222+226
10 | 00000000001111122222222233333333344444444555555555666666777777778888+42
12 | 0011111222333345667779999900001444567888899
14 | 1111133344677889990011222346778
16 | 038900579
18 | 1236825
20 | 235799
22 | 0706
24 | 25
26 | 4

> summary(dis[1,])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  9.066  26.398  35.122  40.466  48.329 264.382 
\end{verbatim}
    Por tanto, IC$_{0.95}$($q$) $\approx$ $[10.5;69.3]$.

    Vamos a recurrir ahora a un procedimiento asintótico.
    Por el teorema central del límite, aproximadamente
    \begin{eqnarray*}
      \bar X
      &\sigue& N\left(\frac12,\frac1{\sqrt{n\cdot(8q+4)}}\right)\\
      {\sqrt{n\cdot(8q+4)}}\cdot(\bar X-1/2)
      &\sigue& N(0,1)\\
      {n\cdot(8q+4)}\cdot(\bar X-1/2)^2
      &\sigue& \chi^2_1\\
      \Pr[{n\cdot(8q+4)}\cdot(\bar X-1/2)^2<\chi^2_{1;1-\alpha}]
      &=&1-\alpha\\
          \Pr\left[{n\cdot(8q+4)}<
          \frac{\chi^2_{1;1-\alpha}}{(\bar X-1/2)^2}\right]
      &=&1-\alpha\\
          \Pr\left[{8q+4}<
          \frac{\chi^2_{1;1-\alpha}}{n\cdot(\bar X-1/2)^2}\right]
      &=&1-\alpha\\
          \Pr\left[{8q}<
          \frac{\chi^2_{1;1-\alpha}}{n\cdot(\bar X-1/2)^2}-4\right]
      &=&1-\alpha\\
          \Pr\left[{q}<
      \frac{\chi^2_{1;1-\alpha}}{8\cdot n\cdot(\bar X-1/2)^2}
      -\frac12\right]
      &=&1-\alpha\\ \ 
          [0;426.86] &\approx& \mbox{IC}_{0.95}(q)
    \end{eqnarray*}
    porque
\begin{verbatim}
> qchisq(.95,1)/(8*n*(mean(x)-1/2)^2) - 1/2
[1] 426.8606
\end{verbatim}
    luego el intervalo sería $[0;426.86]$, demasiado ancho.  La
    presencia de $\bar X-1/2$ en un denominador produce una función
    discontinua que genera valores enormes cuando $\bar X\approx1/2$.
    Sin embargo, la cobertura es buena, por lo que la aproximación
    asintótica parece adecuada; por ejemplo:
\begin{verbatim}
> mean (replicate (100000,
+  33 < qchisq(.95,1)/(8*10*(mean(rbeta(10,33,33))-.5)^2)))
[1] 0.95149
\end{verbatim}
  \end{resoluciOn}
\end{itemize}

\noindent Ten en cuenta que, si $X$ sigue
una distribución $\beta(p,q)$, entonces tiene función de 
densidad $f(x)=x^{p-1}\cdot(1-x)^{q-1}/B(p,q)$, donde $B$ es la
función \texttt{\bfseries beta} en R; esperanza $E(X)=p/(p+q)$ y varianza
$\mathit{Var}(X)=p\cdot q / [(p+q)^2\cdot(p+q+1)]$.  

\paragraph{Problema 2 (3 puntos).}
Se quiere estudiar el tiempo de vida (en años) de los dispositivos eléctricos de la marca Chispazo. Para ello se estudia una muestra de 100 de estos dispositivos eléctricos, recogiendo los tiempos de vida resultantes en el conjunto de datos `Chispazo.RData'. 
\begin{enumerate}
\item[a)] Basándote en los datos recogidos, estudia desde un punto de vista descriptivo si la distribución del tiempo de vida sigue aproximadamente una distribución de Weibull $W(\alpha,\lambda)$ con parámetro de forma $\alpha=2$. Para ello, estima puntualmente el parámetro de escala $\lambda$ de la distribución haciendo uso de su estimador máximo-verosímil.
\item[b)] Suponiendo que el tiempo de vida sigue una distribución de Weibull con parámetro de forma $\alpha=2$, calcula un intervalo de confianza al nivel de confianza $0.99$ para el parámetro de escala $\lambda$ basándote en la distribución asintótica de su estimador máximo-verosímil.
\item[c)] Se quiere realizar un nuevo estudio más preciso que consiga estimar $\lambda$ por medio de un intervalo de confianza asintótico para $\lambda$ al nivel de confianza $0.99$ que tenga una longitud esperada inferior a $0.1$ unidades. Suponiendo que el tiempo de vida sigue una distribución de Weibull con parámetro de forma $\alpha=2$ y haciendo uso de la estimación de $\lambda$ en el estudio preliminar, determina de manera aproximada el tamaño de muestra que se necesita considerar para lograr el objetivo.
\end{enumerate}
\end{document}
