Ver en el Repaso del Práctico 4 el enunciado del teorema de correspondencia entre región de aceptación e intervalos de confianza.
Se da \(X_1, \ldots, X_n\) una poblacion iid con distribucion exponencial de parametro \(\lambda\). Deseamos obtener un intervalo de confianza para la media de la distribucion. Utilizaremos la tecnica de inversion de un test de nivel \(\alpha\) tal que \(H_0 : \lambda = \lambda_0\) vs. \(H_1: \lambda \neq \lambda_0\).
En primer lugar, obtenemos un test de nivel \(\alpha\). Consideremos el test de maxima verosimilitud (LRT) asociado la poblacion dada, que es un test con region de rechazo \(\{ x : \lambda(x) \leq c\}\) para algun \(c\) (a determinar), donde
\[ \Lambda(x) = \dfrac{\sup_{\Theta_0}L(\theta\vert x)}{\sup_{\Theta} L(\theta \vert x)} \] La pdf de la exponencial de parametro \(\lambda\) es \(f(x\vert \lambda) = \lambda^{-1}e^{-x/\lambda}\), \(0 < x < \infty\). La funcion de verosimilitud es \[ L(\theta \vert x_1, \ldots, x_n) = L(\theta \vert x) = f(x\vert \theta) \\ = \prod_{i=1}^n f(x_i \vert \theta) = \lambda^{-n} \prod_{i=1}^n e^{-x_i/\lambda} = \lambda^{-n} e^{-n\bar{x}/\lambda}. \]
Dado el test \(H_0 : \lambda = \lambda_0\) vs. \(H_1: \lambda \neq \lambda_0\), la region del espacio de parametro asociado a \(H_0\), \(\Theta_0\), es el singleton \(\Theta_0 = \{\lambda_0\}\). Por tanto, la optimizacion restringida en \(\Lambda(x)\) (numerador) me da \[ \sup_{\Theta_0}L(\theta\vert x) = \lambda_0^{-n} e^{-n\bar{x}/\lambda_0}. \]
Por otro lado, la optimizacion sin restriccion en \(\Lambda(x)\) (denominador) me da \[ \sup_{\Theta} L(\theta\vert x) = \sup_{\lambda} \lambda^{-n} e^{-n\bar{x}/\lambda}. \] Si \(g(\lambda)\) es la funcion univariada \(g(\lambda) = \lambda^{-n} e^{-n\bar{x}/\lambda}\), entonces \(\log g(\lambda) = -n\log \lambda - n\bar{x}/\lambda\), y \[ g'(\lambda) = -n/\lambda + n\bar{x}/\lambda^2 \Rightarrow g'(\lambda) = 0 \iff \lambda = \bar{x} \] excluyendo el caso \(\lambda = 0\). Verifiquemos que se trata de un maximo: \[ g''(\lambda) = n/\lambda^2 - 2n\bar{x}/\lambda^3 \Rightarrow \\ g''(\lambda=\bar{x}) = n/\bar{x}^2 - 2n/\bar{x}^2 = -n/\bar{x}^2 < 0. \] Por tanto, \[ \sup_{\lambda} \lambda^{-n} e^{-n\bar{x}/\lambda} = \frac{1}{\bar{x}^{n}} e^{-n}. \]
Utlizando los dos resultados anteriores, \[ \Lambda(x) = \dfrac{\lambda_0^{-n} e^{-n\bar{x}/\lambda_0}}{\frac{1}{\bar{x}^{n}} e^{-n}} = \left(\frac{\bar{x}}{\lambda_0}\right)^n e^n e^{-n\bar{x}/\lambda_0}. \]
Para un valor de \(\lambda_0\) fijo, la region de aceptacion es
\[ A(\lambda_0) = \left\{ x : \left(\frac{\sum_i x_i}{n\lambda_0}\right)^n e^n e^{-n\bar{x}/\lambda_0} \geq c\right\} \\ = \left\{x : \left(\frac{\sum_i x_i}{\lambda_0}\right)^n e^{-n\bar{x}/\lambda_0} \geq k_a \right\}. \]
Podemos visualizar esta region de aceptacion en R.
# numero de muestras
n = 8
# representacion de la region de aceptacion
# para cada valor del parametro lambda_0, me doy la region de aceptacion
# s representa la suma de la muestra, sum(x)
A <- function(lambda0){
function(s){(s/lambda0)^n * exp(-s/lambda0)}
}
lambda0 = 3.5 # fijo un valor arbitrario
plot(A(lambda0), 0, 100, xlab=expression(sum(x[i], i==1, n)), main=expression(A(lambda)), ylab="")
# dibujar linear horizontal con el corte en ka
ka = 3000 # valor arbitrario
abline(h=ka)La constante \(k_a\) se debe elegir de modo que
\[ P_{\lambda_0}(X \in A(\lambda_0) = 1-\alpha. \]
Nota. Esta probabilidad se puede calcular haciendo uso de que la suma de v.a. exponenciales se distribuye como una gamma, ie. \(\sum_i X_i \sim \Gamma(n, \lambda)\) y \(\sum_i X_i / \lambda \sim \Gamma(n, 1)\). Esta afirmacion se puede demostrar ya sea directamente utilizando la funcion generatriz de momentos, o por induccion, utilizando que la exponencial es un caso particular de la gamma, ver por ejemplo aqui. Recordar que la pdf de la distribucion \(\Gamma(a, b)\) es \[ f(x) = f(x \vert a, b) = (\Gamma(a)b^a)^{-1} x^{a-1}e^{-x/b} \] para \(0 < x < \infty\) y \(a, b > 0\),
Invertimos la region de aceptacion para obtener el confidence set de tamaño \(1-\alpha\), \[ C(x) = \left\{\lambda : \left(\frac{\sum_i x_i}{\lambda}\right)^n e^{-n\bar{x}/\lambda} \geq k_a \right\}. \]
Visualizamos la region de confianza \(C(x)\) en R.
C <- function(x){
s = sum(x)
function(lambda){(s/lambda)^n * exp(-s/lambda)}
}
# poblacion con distribucion exponencial
# ver help(rexp)
lambda = 5.0
x = rexp(n, rate=1/lambda)
plot(C(x), 0, 40, xlab=expression(lambda), main=expression(C(x)), ylab="")
# dibujar linear horizontal con el corte en ka
ka = 3000 # valor arbitrario
abline(h=ka)A continuación realizamos algunos cálculos analíticos para ilustrar numéricamente el intervalo de confianza obtenido.
La expresion que define a \(C(x)\) depende solamente de \(x\) a traves de la suma de valores, \(\sum_i x_i\). Por tanto, el intervalo de confianza se puede expresar de la manera siguiente:
\[ C\left(\sum_i x_i\right) = \left\{ \lambda: L(\sum_i x_i) \leq \lambda \leq U(\sum_i x_i) \right\}. \] Aqui para determinar a \(L\) y a \(U\), hay dos condiciones:
Si ponemos \(a := \sum_i x_i / L(\sum_i x_i)\) y \(b := \sum_i x_i / U(\sum_i x_i)\), siendo \(a\) y \(b\) constantes, \(a > b\), la condicion 2 queda \[ a^n e^{-a} = b^n e^{-b}. \] Esta ecuacion se puede resolver numericamente (eg. por el metodo de Newton-Raphson con dos variables).
Se trabajara en el caso particular \(n=2\).