\documentclass[12pt]{article}
\usepackage[latin1]{inputenc}
\usepackage[spanish]{babel}
\begin{document}
\title{Ejercicio}
\author{Inferencia}
\maketitle
\raggedright Sea $X$ la variable aleatoria que mide el tiempo, en
meses, entre interrupciones consecutivas de suministro ocurridas en
cierta instalación eléctrica.  Se supone que $X$ sigue una
distribución con densidad $f(x)=\mu^{-1}\cdot e^{-x/\mu}\cdot I(x>0)$
y que las interrupciones son independientes entre sí.  Para estimar la
distribución de $X$ se tomó una muestra aleatoria simple cuyos
resultados fueron los siguientes: \texttt{c(3.35, 1.30, 1.23, 2.13,
  7.70, 11.61, 8.02, 1.44, 3.26, 2.06, 6.06, 0.42, 5.58, 1.40, 6.26,
  0.48, 1.46, 7.33, 7.55, 1.37, 5.41, 0.23, 3.96, 0.86, 0.85, 1.76,
  1.3, 18.36, 2.76, 4.22)}
\begin{enumerate}
\item Calcula el estimador máximo-verosímil de
  $\mu$ y el valor de la estimación para la muestra obtenida.
\item Halla la esperanza y la varianza del estimador máximo-verosímil de 
  $\mu$.
\item Obtén un intervalo de confianza con coeficiente $1-\alpha=0.95$
  para la media y la mediana de $X$.
% \item Contrasta al nivel de significación $\alpha=0.01$, las hipótesis
% $$H_0: X \equiv \; Exponencial \qquad H_1:X \ne \; Exponencial$$
% \item Contrasta al nivel de significación $\alpha=0.01$, las hipótesis
% $$H0: \mu \le 3 \qquad H_1:\mu > 3$$
\item Calcula la estimación máximo-verosímil de $\mu$, suponiendo que
  los datos de la muestra estuvieran agrupados en las siguientes
  clases:
%$$[0,\, 1),\, [1,\, 1.5),\, [1.5, \,3),\, [3, 4.5),\, [4.5,\,  7.5), (\ge 7.5) $$

 \begin{center}
 \begin{tabular}{c|c|c|c|c|c|c}
 Intervalo & $[0;1)$ & $[1;1.5)$ & $[1.5;3)$ & $[3;4.5)$ & $[4.5;7.5)$ & $\ge 7.5$\\
 \hline Frecuencia & $5$ & $7$ & $4$ & $4$ & $5$ & $5$
 \end{tabular}
 \end{center}



% \item Contrasta al nivel de significación $\alpha=0.01$, las hipótesis
% $$H_0: X \equiv \; Exponencial \qquad H_1:X \ne \; Exponencial$$ con la muestra agrupada en las clases del apartado anterior.

% \end{enumerate}

% Solución. \\

% \begin{enumerate}
% \item La función de densidad de la variable es 
% $$ f(x) = \lambda \, e^{-\lambda \, x} I(x> 0) = 
% \frac{1}{\mu} \, e^{- x/ \mu} I(x> 0)$$
% donde $\lambda$ representa la \textit{velocidad} con la que ocurren los cortes de suministro y $\mu$ la media del tiempo transcurrido entre dos consecutivos, con $E(X) = \mu= \frac{1}{\lambda}$.

% La función de verosimilitud es:
% $$L(x_1, \ldots, x_n, \mu)= \prod_{i=1}^n \, \frac{1}{\mu} \, e^{- x_i/ \mu} I(xi> 0) = \frac{1}{\mu^n } \, e^{- \sum_{i=1}^n x_i/ \mu} \prod_{i=1}^n  I(xi> 0)$$

% Utilizando la transformación $ln \left(L(x_1, \ldots, xn, \mu) \right)$, derivando respecto a $\mu$ e igualando a cero, se facilita el cálculo del valor donde se alcanza el máximo de la verosimilitud.

% %El estimador máximo verosímil es el valor de $\mu$ que hace máxima la expresión anterior, que coincide con 
\end{enumerate}

% $$ln \left(L(x_1, \ldots, xn, \mu) \right) = - n\, ln(\mu) -  \frac{ \sum_{i=1}^n x_i}{\mu}$$ 
% $$ \frac{d\, ln \left(L(x_1, \ldots, xn, \mu) \right)}{d \, \mu}= -\frac{n}{\mu} + \frac{ \sum_{i=1}^n x_i}{\mu^2}$$ 

% Igualando a cero se obtiene:
% $$0= \frac{n}{\hat{\mu}} + \frac{ \sum_{i=1}^n x_i}{\hat{\mu}^2}
% \Longleftrightarrow  n \, \hat{\mu}=  \sum_{i=1}^n x_i$$
% El estimador máximo verosímil de la media teórica es la media muestral: 
% $$\hat{\mu} = \frac{\sum_{i=1}^n x_i}{n} $$ 

% \begin{verbatim}
%  ## Librerías
% import numpy as np
% from scipy import stats as st
% import pandas as pd

% ## muestra aleatoria simple
% x= [3.35, 1.30, 1.23, 2.13, 7.70, 11.61, 8.02, 1.44, 3.26, 2.06, 6.06, 
%     0.42, 5.58, 1.40, 6.26, 0.48, 1.46, 7.33,  7.55, 1.37, 5.41, 0.23,
%     3.96, 0.86, 0.85, 1.76,1.30, 18.36, 2.76, 4.22]
% n= len(x); print("tamaño de la muestra=", n) ## tamaño de la muestra
% ## a) Estimación máximo verosímil
% mm=np.mean(x)## media muestral
% print("Estimador máximo verosímil, mm.round(4))

% ## c) Intervalo de confianza exacto
% alfa= 0.05 ## 1-alfa = coeficiente de confianza
% n= len(x); print(n) ## tamaño de la muestra
% \end{verbatim}

% \item Calcula la esperanza y la varianza del estimador.\\

% $$E(\hat{\mu}) = E\left(\frac{\sum_{i=1}^n x_i}{n}\right) = 
% \frac{\sum_{i=1}^n E(X_i) }{n}= \frac{\sum_{i=1}^n \mu }{n}= \mu$$

% $$Var(\hat{\mu}) = Var\left(\frac{\sum_{i=1}^n x_i}{n}\right) = 
% \frac{Var\left(\sum_{i=1}^n X_i\right) }{n^2}=
% \frac{\sum_{i=1}^n Var(X_i) }{n^2}= \frac{\sum_{i=1}^n \sigma^2 }{n^2}= \frac{\sigma^2}{n}= \frac{\mu^2}{n}$$

% Ya que por ser las variables $X_i$ independientes $Var\left(\sum_{i=1}^n X_i\right) = \sum_{i=1}^n Var(X_i) $. Además por seguir una distribución exponencial $Var(X_i)= \mu^2= \frac{1}{\lambda^2}$

% \item Obtén un intervalo de confianza con coeficiente $1-\alpha=0.95$ para la media y la mediana de $X$.

% El intervalo de confianza se va a calcular utilizando el método de la función pivote y el método bootstrap.

% \begin{itemize}
%  \item La función pivote se construye usando la media muestral y utilizando la siguiente propiedad de la distribucción exponecial:
% $$ \mbox{Si } X\equiv Exp(\mu) \Rightarrow Y=\frac{X}{\mu} \equiv Exp(1)$$
% Los valores que puede tomar $Y$ son los mismos que X, es decir el intervalo $(0, \infty)$. Además, la función de distribución de $Y$ es:
% $$F(y) = P(Y\le y)  = P(\frac{X}{\mu} \le y) = P(X \le \mu \, y) = 
% \int_0 ^{\mu \, y} \, f(x) \, dx = $$
% $$\int_0 ^{\mu \, x} \,\frac{1}{\mu}  \, e^{ - x/ \mu} dx = 1- e^{ - \mu \, y \, / \mu} = 1- e^{-y}  $$

% Por tanto la función de densidad vale:
% $$f(y) = \frac{d \, F(y)}{y} = e^{-y} I(y > 0) \equiv Exp(1) $$

% En consecuencia $$\frac{\bar{X}}{\mu}= \frac{1}{n} \sum_{i=1}^{n} \frac{X_i}{\mu} = \frac{1}{n} \sum_{i=1}^{n} Yi \equiv \frac{1}{n} \sum_{i=1}^{n} Exp(1) =  \frac{1}{n} Gamma(n,1) = Gamma(n,\frac{1}{n}) $$

% El intervalo con coeficiente de confianza 0.95 se calcula de la siguiente manera:

% $$0.95 = P(a \le \frac{\bar{X}}{\mu} \le b) = 
% P(\frac{1}{b} \le \frac{\mu}{\bar{X}} \le \frac{1}{a}) = 
% P(\frac{\bar{X}}{b} \le \mu \le \frac{\bar{X}}{a})$$ 

% Para obtener los valores de $a$ y $b$ se plantean las siguientes ecuaciones:

% $$ P\left(Gamma(n, \frac{1}{n}) \le a \right) = \frac{0.05}{2}= 0.025;  \qquad 
%    P(Gamma(n, \frac{1}{n}) \le b) = 1- \frac{0.05}{2}= 0.975$$
   
%    El intervalo de confianza exacto de la media es:
   
%   \begin{itemize}
%      \item IC exacto (media)= [2.8745; 5.9148]
%      \end{itemize}
% %   $$IC(media)= [2.8745, 5.9148]$$

%    \item Intervalo de confianza bootstrap básico\\
%    Los valores de $a$ y $b$ se estiman a partir de las réplicas bootstrap del estadístico $\bar{X}^b - \bar{X}$, para que se verifique
%    $P(a \le \bar{X}^b - \bar{X} \le b)= 1-\alpha$
   
%    \vspace{0.3cm}
   
%    A continuación se muestran los intervalos de confianza bootstrap por el método básico y percentil con la librería \textit{bootstrap de scipy.stats}. 
%    \begin{itemize}
%      \item IC boostrap básico (python)= ( 2.4772 ;  5.2524 )
%      \item IC boostrap percentil** (python)=( 2.7126 ;  5.5168 )
%      \item intervalo bootstrap básico (programado)=  [2.4513; 5.2504]
%      \item intervalo bootstrap percentil** (programado)=  [2.731; 5.53]
% \end{itemize}

%    \end{itemize}

% \vspace{0.3cm}   
   
% \begin{verbatim}
% ## Librerías
% ## c) Intervalo de confianza exacto
% alfa= 0.05 ## 1-alfa = coeficiente de confianza
% n= len(x); print(n) ## tamaño de la muestra

% ## c.1) Función pivote. IC= (mm/b, mm/a)
% a= st.gamma.ppf(alfa/2, n,loc=0, scale=1/n) 
% b= st.gamma.ppf(1-alfa/2, n,loc=0, scale=1/n)
% ei= mm/b; es= mm/a
% print("IC(media)=", [ei.round(4), es.round(4)])

% ## c.1) Intervalos de confianza bootstrap 
% ## Con la librería bootstrap de scipy.stats
% from scipy.stats import bootstrap
% xb= x ## almacena las muestras boostrap
% xb = (xb,)  # samples must be in a sequence
% nrb= 9999 ## número de réplicas bootstrap
% ## boostrap básico
% res_basic = bootstrap(xb, np.mean, confidence_level=1-alfa,
%             n_resamples=nrb, method='basic', random_state=rng)
% ## boostrap percentil
% res_perc = bootstrap(xb, np.mean, confidence_level=1-alfa,
%            n_resamples=nrb, method='percentile', random_state=rng)
% ## intervalos
% ei_bb= res_basic.confidence_interval.low
% es_bb= res_basic.confidence_interval.high
% ei_bp= res_perc.confidence_interval.low
% es_bp= res_perc.confidence_interval.high
% print("IC boostrap básico= (",ei_bb.round(4),"; ",es_bb.round(4),")")
% print("IC boostrap percentil=(",ei_bp.round(4),"; ",es_bp.round(4),")")

% ## itervalos bootstrap programado.
% mm=np.mean(x) ## media muestral
% b_per=[]
% dif=[]
% for j in range(1, nrb+1):
%     xb=np.random.choice(data[0], n) # muestra bootstrap
%     mmb=np.mean(xb)## media muestral bootstrap
%     dif.append(mmb-mm)
%     b_per.append(mmb)

% ## intervalo bootstrap básico
% cuan= np.quantile(dif, [alfa/2, 1-alfa/2], axis=0)
% eib= mm - cuan[1] ; esb= mm - cuan[0] 
% print("intervalo bootstrap básico= ",[eib.round(4), esb.round(4)])

% ## intervalo bootstrap percentil
% cuan2= np.quantile(b_per, [alfa/2, 1-alfa/2], axis=0)
% eib2= cuan2[0] ; esb2= cuan2[1] 
% print("intervalo bootstrap percentil= ",[eib2.round(4),esb2.round(4)])
% \end{verbatim}

% \item Contrasta al nivel de significación $\alpha=0.01$, las hipótesis
% $$H0: X \equiv \; Exponencial \qquad H_1:X \ne \; Exponencial$$
% La distribución exponencial tiene un parámetro desconocido, que es la media teórica, $\mu$, y por tanto se debe de emplear el test de Kolmogorov-Smirnov, con la corrección de Lilliefors

% El resultado del test es: 
% \begin{center}
%  Lilliefors:  estadistico= 0.0986 ; pvalor= 0.7975
% \end{center}
% Código de python:
% \begin{verbatim}
%  ## Test de bondad de ajuste de Lilliefors
%  import statsmodels.stats.diagnostic as ssdiag ## Test de Lilliefors
% ## H0: X= Exponencial; H1: X= no sigue una Exponencial
% lillie=ssdiag.lilliefors(x, dist='exp', pvalmethod='table')
% print("Lilliefors: ","estadistico=",round(lillie[0],4), 
%       "; pvalor=", round(lillie[1],4))  
% \end{verbatim}


% \item Contrasta al nivel de significación $\alpha=0.01$, las hipótesis
% $$H0: \mu= 3 \qquad H_1:\mu > 3$$

% Bajo $H_0$ las muestras deberían tener una media muestral más pequeña que si fuera cierta la hipótesis $H_1$. Por tanto una región crítica \textit{razonable} sería la siguiente:
% $$RC= \{(x_1, \ldots, x_n) | \bar{x} > k \}$$
% Dado que la media de la muestra obtenida vale $\bar{x}= 3.9907$ el p-valor del contraste es:
% $$P(  \bar{X} > 3.9907 |H_0) = 
% P(\frac{1}{n} \, \sum_{i=1}^{n} X_i > 3.9907 |H_0) = 
% P(\frac{1}{n} \, \sum_{i=1}^{n} Exp(3) > 3.9907) =$$
% $$P(\frac{1}{n} \,Gamma(n, 3) > 3.9907)=
% P(Gamma(n, \frac{3}{n}) > 3.9907) = 0.0445$$

% Al nivel de significación $\alpha=0.05$ se rechazaría la hipótesis nula, pero dado que el nivel de significación elegido es 0.01 la decisión que se debe tomar es no rechazar que la media teórica valga tres.

% \begin{verbatim}
% ## X= Exponencial. Contrastar H0: media=3; H1: media > 3
% ## pvalor= P(Gamma(n,loc=0, scale=3/n) > media muestral)
% pvalor= 1- st.gamma.cdf(mm, n,loc=0, scale=3/n) 
% print("Contraste: H0: media=3; H1: media > 3. pvalor=", pvalor.round(4))
% \end{verbatim}

% \item Calcula la estimación máximo-verosímil de  $\mu$ con los datos  agrupados en clases.

% La probabilidad del intervalo $C_j$, cuyos extremos son $x_{j-1}$ y $x_j$ vale:
% $$P(C_j) = \int_{x_{j-1}}^{x_j} \, f(x) \, dx =
% \int_{x_{j-1}}^{x_j} \, \frac{1}{\mu} e^{-x/ \mu} \, dx =
% e^{-x_{j-1}/ \mu} - e^{-x_{j}/ \mu} $$

% En consecuencia la función de verosimilitud es:
% $$L(C_1, \ldots, C_n, \mu) = \prod_{i=1}^n P(C_i) = 
% \prod_{j=1}^k P(C_j)^{n_j} = 
% \prod_{j=1}^k e^{-x_{j-1}/ \mu} - e^{-x_{j}/ \mu} $$

% $$log(L(C_1, \ldots, C_n, \mu)) = 
% \sum_{j=1}^k n_j log(e^{-x_{j-1}/ \mu} - e^{-x_{j}/ \mu}) $$

% donde $k= 6$, $n_j=$ las frecuencias de cada clase $(j=1,..., k)$ y los extremos de los intervalos valen 
% $x_0=0, x_1=1, x_2=1.5, x_3=3; x_4=4.5, X_5= 7.5, x_6= \infty$
% En esta función no se puede obtener el máximo de manera analítica y se debe resolver por métodos numéricos.

% $$EMV (\mbox{estimación con datos agrupados})=  4.053$$

% \begin{verbatim}
% ## Estimación máximo verosímil con los datos agrupados
% O= np.array([5, 7, 4, 4, 5, 5]) ## frecuencias observadas
% ## P(Cj) = P(X <= Xj) - F(X <= Xj-1)
% cortes= [0, 1, 1.5, 3, 4.5, 7.5] ## extremos inferiores de los intervalos
% ## Función a minimizar: - función de verosimilitud
% def f(l):
%     p2= np.array(list(st.expon.cdf(cortes[1:6], loc=0, scale=l)) + [1])
%     p1= np.array(st.expon.cdf(cortes, loc=0, scale=l))
%     prob= p2-p1
%     return - np.sum(O*np.log(prob))
% f(2)  ## prueba de la función
% ## Gráfica de - log(Verosimilitud)
% para = np.linspace(0.5, 12, 1000 )
% y=[]
% for v in para:
%     y.append(f(v))
% ## gráfica    
% pl.plot(para, y)
% pl.xlabel("media")
% pl.ylabel("- log(Verosimilitud)")
% indice= np.array(y).argmin(); 
% print("aproximación del máximo verosímil=", para[indice].round(4))
% ## estimación numérica del máximo verosímil
% inicio= (cortes[3]+cortes[4])/2; 
% print("Valor de inicio del mínimo= ",inicio) 
% salida= optimize.minimize(f, x0= inicio, bounds=[(0.5, 12)]); 
% emv_a = salida.x[0] ## estimación máximo verímil con datos agrupados
% print("EMV numérico (datos agrupados)= ", emv_a)
% \end{verbatim}

% \item Contrasta al nivel de significación $\alpha=0.01$, las hipótesis $$H0: X \equiv \; Exponencial \qquad H_1:X \ne \; Exponencial$$

% Como la muestra está agrupada en intervalos se debe utilizar el test chi cuadrado de bondad de ajuste, cuyo estadístico viene dado por:

% $$ W= \sum_{j=1}^k \frac{(O_j - E_j)^2}{E_j}$$
% donde $O_j$ son las frecuencia observadas de cada intervalo y $E_j= n \ P(C_j)$ las frecuencias esperadas bajo$H_0$

% En este caso, la tabla de frecuencias observada, probabilidades y frecuencias esperadas es la siguiente:
% \begin{center}
% \begin{tabular}{c|c|c|c|c|c|c}
% Intervalo & [0,\,1) & [1,\,1.5) & [1.5,\,3) & [3,\, 4.5) & [4.5,\, 7.5) & $\ge 7.5$\\
% \hline
% O & 5 & 7 & 4 & 4 & 5 & 5\\
% \hline
% P & 0.2186 & 0.0907 & 0.2136 & 0.1476 & 0.1723 & 0.1572\\
% \hline
% E & 6.5580 & 2.7220& 6.4094 & 4.4267 & 5.1691 & 4.7148
% \end{tabular}
% \end{center}

% Al realizar el contraste de hipótesis, se obtiene:
% \begin{center}
%  Estadístico=  8.0723 ; grados libertad= 4 ; pvalor= 0.911
% \end{center}

% Aunque las frecuencia esperada en el intervalo $[1, \, 1.5)$ es un poco baja, está cerca de tres y la aproximación a una distribución chi-cuadrado es aceptable. Por otra parte, el pvalor del contraste indica que es un valor muy frecuente bajo $H_0$ y no hay ninguna evidencia para  rechazar la hipótesis nula de que $X$ siga una distribución exponencial

% \begin{verbatim}
% ## Contrastar la hipótesis: H0: X= Exponencial; H1: X= no Exponencial
% ## Test chi- cuadrado de bondad de ajuste.
% p2= np.array(list(st.expon.cdf(cortes[1:6], loc=0, scale=emv_a)) + [1])
% p1= np.array(st.expon.cdf(cortes, loc=0, scale=emv_a))
% prob= p2-p1 ## probabilidades de los intervalos
% E= np.array(n*prob) ## frecuencias esperadas bajo H0
% print("frecuencias esperadas", E)
% W= sum((O-E)**2/E) ## estadístico chi-cuadrado
% ## grados de libertad = len(O) - 1 - número de parámetros estimados
% pvalor= st.chi2.cdf(W, len(O)-2) ## pvalor del contrate
% print("Estadístico= ", W.round(4), "; grados libertad=", len(O)-2, 
%       "; pvalor=", pvalor.round(4))
% \end{verbatim}
\end{document}
