Práctica 1: Inferencias de Wald y basadas en la RV

Inferencias de Wald: IC y CH

Sean \(X_1,\dots, X_n\) v.a.i.i.d. \(P_\theta, \theta\in \Theta\subset\mathbb{R}; f(x,\theta).\)

Como hemos visto, en situaciones regulares el EMV (\(\hat{\theta}\)) para \(\theta\) es CAN y AE: \[\begin{equation} \sqrt{n}(\theta-\hat{\theta})\xrightarrow{\mathcal{L}} N(0,V^2(\theta)), \quad (Eq. 1) \end{equation}\] donde \(V^2(\theta)=\frac{1}{I_1(\theta)}.\)

Las inferencias acerca de \(\theta\) basadas en esta distribución asintótica se conocen como inferencias de Wald. Por un lado, para obtener un IC para \(\theta,\) se necesita estimar \(I_1(\theta).\) Tal y como vimos la IF observada, que definimos como \(\widehat{I_1(\theta)}=-\frac{\partial^2}{\partial \theta^2}\log f(X,\theta)\vert_{\theta=\hat{\theta}}\), nos permite estimar esta cantidad, y concretamente: \[\widehat{V^2(\theta)}=\hat{V}^2(\theta)=\frac{1}{\widehat{I_1(\theta)}},\]

es decir, \(V^2(\theta)\) se estima como la inversa de la IF observada, y de donde también se sigue que:

\[\sqrt{n}(\theta-\hat{\theta})\xrightarrow{\mathcal{L}} N(0,\widehat{V^2(\theta)})\text{, o equivalentemente, } \hat\theta\cong N(0,\frac{\widehat{V^2(\theta)}}{n})\] NOTA: Cuidado con estas dos expresiones, soleis cometer fallos con la \(n.\)

Por tanto, \(Z_0=\frac{\sqrt n|(\hat{\theta}-\theta)|}{\sqrt{\hat{V}^2(\theta)}}\cong N(0,1).\) Así, fijado \(\alpha\), un IC un Wald \((1-\alpha)\%\) para \(\theta\) queda definido por: \[P_{\theta}( \frac{\sqrt n|(\hat{\theta}-\theta)|}{\sqrt{\hat{V}^2(\theta)}}\leq\lambda_{1-\alpha/2})\simeq 1-\alpha,\] donde \(\lambda_{1-\alpha/2}\) es el cuantil \(1-\alpha/2\) de una normal estandar (\(qnorm(1-\alpha/2)\) en {R}). O equivalentemente, despejando \(\theta\) en la expresión anterior se obtiene el IC de Wald pedido: \[P_{\theta}(\theta\in[\hat{\theta}\pm qnorm(1-\alpha/2)\frac{\sqrt{\hat{V}^2(\theta)}}{\sqrt{n}}])\simeq 1-\alpha\] Gráficamente:

El resultado anterior de la Eq. 1 también es clave para resolver CH para \(\theta\): \[\begin{equation} H_0:\theta=\theta_0 \quad \textit{vs}\quad H_1:\theta\neq \theta_0. \end{equation}\].

Tal y como vimos el estadístico de Wald: \[Q_W=\frac{n(\hat{\theta}-\theta_0)^2}{\widehat{V^2(\theta)}}\simeq n(\hat{\theta}-\theta_0)^2 \widehat{I_1(\theta)}\overset{H_0}{\cong}\chi^2_{1}\] por tratarse del cuadrado de una distribución normal estándar. Esta distribución nos permitirá calcular el p-valor como \(P_{\theta_0}(Q_W>q_{W,obs}).\)

Inferencias basadas en el estadístico RV: IC y CH

Sean \(X_1,\dots, X_n\) v.a.i.i.d. \(P_\theta, \theta\in \Theta\subset\mathbb{R}; f(x,\theta).\) Dado un nivel de significación \(\alpha\) se quiere resolver el CH para \(\theta\): \[\begin{equation} H_0:\theta=\theta_0 \quad \textit{vs}\quad H_1:\theta\neq \theta_0. \end{equation}\] Según lo estudiado, el estadístico de la RV se define como: \[Q_L(\boldsymbol X)=2[\log L(\hat{\theta},\boldsymbol X)-\log L(\theta_0,\boldsymbol X)]\] y es asintóticamente equivalente a \(Q_w.\) Luego el test de RV de nivel \(\alpha\) para resolver CH para \(\theta\) tiene una región crítica de tamaño \(\alpha\) aproximadamente que viene dada por \((Q_L(X)>\chi^2_{1,1-\alpha})\simeq (Q_L(X)>qchisq(1,1-\alpha)),\) donde \(\chi^2_{1-\alpha,1}\) es el cuantil \(1-\alpha\) de una \(\chi^2_1.\) El p-valor del test vendrá dado por \(P_{\theta_0}(Q_L>q_{L,obs}).\)

De la relación entre test de nivel \(\alpha\) e IC \((1-\alpha)\), se tiene que el IC basado en el estadístico de la RV (ICRV\((1-\alpha)\)) se corresponde con el conjunto de valores de \(\theta\) para los que el test de la RV de nivel \(\alpha\) no rechaza la \(H_0\) es decir: \[\begin{align} ICRV(1-\alpha)&=\lbrace \theta_0: Q_L(X) \leq qchisq(1-\alpha,1) \rbrace\\ &= \lbrace \theta_0: 2[\log L(\hat{\theta},X)-\log L(\theta_0,X)] \leq qchisq(1-\alpha,1) \rbrace\\ &= \lbrace \theta_0: \log L(\theta_0,X) \geq \log L(\hat{\theta},X)-\frac{qchisq(1-\alpha,1)}{2} \rbrace\\ &= \lbrace \theta_0: \log L(\theta_0,X) \geq d_1 \rbrace \end{align}\]

Ejercicio 1: Modelo de Poisson.

Considere \(X_1,\dots,X_{50}\) m.a.s. de \(\mathcal{P}(\lambda),\) \(f(x,\lambda)=e^{-\lambda}\frac{\lambda^x}{x!},x=0,1,2,\dots.\) Para los datos observados se tiene \(S=\sum_{i=1}^{50} P_{\lambda}(X_i=x_i)=75.\)

1. Obtener y representar la función de verosimilitud y log verosimilitud.

Verosimilitud: \(L(\lambda,X_1,\dots,X_n)\propto e^{-n\lambda}\frac{\lambda^{\sum_{i=1}^n X_i}}{\prod_{i=1}^n X_i!}\)

Log verosimilitud: \(\log L(\lambda,X_1,\dots,X_n) = -n\lambda +\sum_{i=1}^n X_i \log \lambda\)

#Datos
n<-50
S<-75
#Verosimilitud X1,...x_n ~ P(lambda)
#Todo aquello que no depende de lambda paodemos suprimirlo
v<-function(lambda){
  return(exp(-n*lambda)*lambda^S)
}
#Log-verosimilitud
lv<-function(lambda){
  return(-n*lambda+S*log(lambda))
}
par(mfrow=c(1,2))
plot(seq(0,10,length.out=500),v(seq(0,10,length.out=500)),xlab="lambda",ylab="L(lambda)",
     main="Verosimilitud. n=50, S=75",type="l")
plot(seq(0,10,length.out=500),lv(seq(0,10,length.out=500)),xlab="lambda",ylab="l(lambda)",
     main="Log verosimilitud. n=50, S=75",type="l")

2. Calcular el EMV de forma analítica. Ecuación de verosimilitud: \(\frac{\partial}{\partial \lambda}\log L(\lambda,X_1,\dots,X_n)=-n+\frac{\sum_{i=1}^n X_i}{\lambda}=0,\) de donde se sigue \(\hat{\lambda}=\overline{X}=\frac{\sum_{i=1}^n X_i}{n}=\frac{75}{50}=1.5\)

#EMV
emv<-S/n
emv

## [1] 1.5

par(mfrow=c(1,2))
plot(seq(0,10,length.out=500),v(seq(0,10,length.out=500)),xlab="lambda",ylab="L(lambda)",
     main="Verosimilitud. n=50, S=75",type="l")
abline(v=emv,col=2)
plot(seq(0,10,length.out=500),lv(seq(0,10,length.out=500)),xlab="lambda",ylab="l(lambda)",
     main="Log verosimilitud. n=50, S=75",type="l")
abline(v=emv,col=2)

3. Obetener el EMV de forma numérica (aproximada) mediante la función {optim}. No siempre es posible calcular EMV de forma analítica y es común recurrir a métodos numéricos (por ejemplo Newton Raphson). En {R} podemos usar la función {optim} que minimiza funciones lineales y no lineales usando Nelder–Mead, quasi-Newton o de descenso de gradiente.

#Leer https://cran.r-project.org/doc/contrib/Optimizacion_Matematica_con_R_Volumen_I.pdf

#optim(par, fn, gr = NULL, ...,
#      method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN",
#                 "Brent"),
#      lower = -Inf, upper = Inf,
#      control = list(), hessian = FALSE)
#By default optim performs minimization, but it will maximize negative functions.


#par: Initial values for the parameters to be optimized over.

#fn: A function to be minimized, with first argument the vector of parameters over which minimization is to take place. It should return a scalar result.

#method: The default method is an implementation of that of Nelder and Mead. Method "Brent" is for one-dimensional problems only

#lower, upper: Bounds bounds in which to search for method "Brent".

#hessian: Logical. Should a numerically differentiated Hessian matrix be returned? Valor de la segunda derivada de fn
  
#Caso uniparamétrico basta con:
mlv<-function(lambda){
  return(-lv(lambda))
}
minimo<-optim(par=3, fn=mlv, method = "Brent",
      lower = 0.1, upper = 3,  hessian = TRUE)

#> minimo$par
#> [1] 1.5
#> 
#> minimo$value
#> [1] 44.59012
#> 
#> minimo$hessian   
#> [1] 33.33336
par(mfrow=c(1,1))
plot(seq(0,10,length.out=500),mlv(seq(0,10,length.out=500)),xlab="lambda",ylab="l(lambda)",
     main="-Log verosimilitud. n=50, S=75",type="l")
abline(v=emv,col=2)

emv<- minimo$par
emv

## [1] 1.5

NOTA: Estamos trabajando con {mlv}, de forma que la inversa de {minimo$hessian} es un estimador de la varianza asintótica.

4. Calcular IC de Wald para \(\lambda\) con confianza aproximada del \(95\%.\)

Podemos llegar de dos maneras diferentes.

Distribución asintótica del EMV (TCL): \(\sqrt{n}(\overline{X}-\lambda)\xrightarrow{\mathcal{L}}N(0,\lambda)\), puesto que \(E_\lambda(X_i)=\lambda\) y \(Var_\lambda(X_i)=\lambda \forall i.\)

Distribución asintótica del EMV (Eq. 1): \(\sqrt{n}(\overline{X}-\lambda)\xrightarrow{\mathcal{L}}N(0,\frac{1}{I_1(\lambda)}),\) donde \(I_1(\lambda)=-E_\lambda(\frac{\partial^2}{\partial \lambda^2}\log f(X,\lambda))=-E_\lambda(\frac{-X}{\lambda^2 })=\frac{1}{\lambda}.\)

Para poder hacer inferencias, necesitamos estimar \(I_n(\lambda)\) mediante la IF observada: Calculamos \(\frac{\partial^2}{\partial\lambda^2}\log L(\lambda,X_1,\dots,X_n)=-\frac{\sum_{i=1}^n X_i}{\lambda^2}\), \(\widehat{I_n(\lambda)}=-\frac{\partial^2}{\partial \lambda}\log f(\boldsymbol{X},\lambda)\vert_{\theta=\hat{\lambda}}=\frac{\sum_{i=1}^n X_i}{\lambda^2}\vert_{\lambda=\overline{X}}=\frac{n^2}{\sum_{i=1}^n X_i}=\frac{2500}{75}.\)

Y también se tiene: \((\overline{X}-\lambda)\xrightarrow{\mathcal{L}}N(0,\frac{75}{2500}),\) donde \(se(\hat{\lambda})=\sqrt{\frac{75}{2500}}.\)

Por tanto, el IC \((1-\alpha)\%\) de Wald para \(\lambda\) es: \[(\overline{X}\pm qnorm(0.975)\sqrt{\frac{75}{2500}})\simeq(1.5\pm qnorm(0.975)\sqrt{0.03})\simeq(1.16,1.84)\]

#var asintotica del emv
emv.var<-1/minimo$hessian # la inversa de la salida Hessian
emv.var

##            [,1]
## [1,] 0.02999997

#IC Wald con confianza 0.95
emv-qnorm(0.975)*sqrt(emv.var)

##          [,1]
## [1,] 1.160524

emv+qnorm(0.975)*sqrt(emv.var)

##          [,1]
## [1,] 1.839476

5. Calcular el estadístico de Wald para contrastar \(H_0:\lambda=1 \quad \textit{vs}\quad H_1:\lambda\neq 1\). ¿Cuál es el p-valor del test?

\(Q_W=\frac{n(\hat{\lambda}-\lambda_0)^2}{\widehat{V^2(\lambda)}}= n(\hat{\lambda}-\lambda_0)^2 \widehat{I_1(\lambda)}=\frac{n(\overline{X}-1)^2}{\overline{X}}\overset{H_0}{\cong}\chi^2_{1}\)

Región crítica: \((Q_W>C_{0.05}),\) donde \(C_{0.05}=qchisq(0.95,1)=3.841,\) pues \(Q_W\sim_{H_0} \chi^2_1.\)

P-valor: \(P_{\lambda=1}(Q_W>q_{W,obs})\cong 1-P_{\lambda=1}(\chi^2_1\leq q_{W,obs})=1-pchisq(q_{W,obs},1),\) donde \(q_{W,obs}=\frac{n(\overline{X}-1)^2}{\overline{X}}=\frac{50(1.5-1)}{1.5^2}\)

#Estadistico Wald al tipicar la distribuciaon asintótica
Q_Wobs<-((emv-1)^2)/emv.var
Q_Wobs

##          [,1]
## [1,] 8.333341

#Region crítica
qchisq(0.95,1)

## [1] 3.841459

#P-valor
pvalorW<-1-pchisq(Q_Wobs,1)
pvalorW # Se rechaza H_0

##             [,1]
## [1,] 0.003892401

6. Calcular el estadístico test de la RV para contrastar \(H_0:\lambda=1 \quad \textit{vs}\quad H_1:\lambda\neq 1\). ¿Cuál es el p-valor del test? Estadístico RV: \[\begin{align} Q_L&=2[\log L(\hat\lambda,X_1,\dots,X_n)-\log L (\lambda_0,X_1,\dots,X_n)]\\ &=2[\log L(\overline{X},X_1,\dots,X_n)-\log L(1,X_1,\dots,X_n)]\\ &=2[-n\overline{X}+n\overline{X}\log\overline{X}+n]\overset{H_0}{\cong}\chi^2_{1} \end{align}\]

P-valor: \(P_{\lambda=1}(Q_L>q_{L,obs})\cong 1-P_{\lambda=1}(\chi^2_1\leq q_{L,obs})=1-pchisq(q_{L,obs},1),\) donde \(q_{L,obs}=2[-75+75\log(1.5)+50]\)

#Estadistico RV
Q_Lobs<-2*(-minimo$value+n)#donde toma el valor minimo
Q_Lobs#DISTINTO DE Q_W

## [1] 10.81977

#P-valor
pvalorL<-1-pchisq(Q_Lobs,1)
pvalorL # Se rechaza H_0, si est Wald ha rechazado, este también

## [1] 0.001004222

7. Calcular IC basado en la RV para \(\lambda\) con confianza aproximada del \(95\%.\) \[\begin{align} ICRV(0.95)&=\lbrace \lambda_0: 2[\log L (\overline{X},\boldsymbol{X})-\log L (\lambda_0,\boldsymbol{X})]\leq qchisq(0.95,1) \rbrace\\ &=\lbrace \lambda_0: \log L(\lambda_0,\boldsymbol{X})\geq\log L(\overline{X},\boldsymbol{X})-\frac{qchisq(0.95,1)}{2}\rbrace\\ &=\lbrace \lambda_0: \log L(\lambda_0,\boldsymbol{X})\geq d_1\rbrace \end{align}\]

#Calcular la constante d_1 
d1<--minimo$value-qchisq(0.95,1)/2
d1

## [1] -46.51085

#Representacion de log verosimilitud en un entorno de emv, p.e. 0.5 y 3
bound<-seq(0.5,3,length.out=100)
plot(bound,lv(bound),type="l")
#añado la constante en d
abline(h=d1)

#ahora queda buscar esos puntos 
#n un par de intentos se puede encontrar 
#abline(v=1.1)
#abline(v=1.2)
abline(v=1.19)

#abline(v=1.8)
#abline(v=1.9)
abline(v=1.86)

#[1.86,1.19]

8. Calcular IC de Wald para \(g(\lambda)=P_\lambda(X=0)=e^{-\lambda}\) con confianza aproximada del \(95\%.\) Por un lado, el EMV es invariante para cualquier función \(g\) (con inversa), se tiene que \(\widehat{g(\lambda)}=g(\hat{\lambda})=g(\overline{X})=e^{-\overline{X}}.\) Por otro lado, el Delta Método permite a partir de la distribución asintótica del EMV calcular la distribución de cualquier función \(g\) (derivable) del EMV como: \[\sqrt{n}(g(\hat\lambda)-g(\lambda))\xrightarrow{\mathcal{L}}N(0,\widehat{V^2(\lambda)}(g^{'}(\hat\lambda))^2,\] donde \(g(\hat\lambda)=e^{-\overline{X}},\) \(g(\lambda)=e^{-\lambda},\) \(\widehat{V^2(\lambda)}=\overline{X}\) y \((g^{'}(\hat\lambda))^2=(-e^{-\hat\lambda})^2=e^{-2\overline{X}}.\) De donde se sigue: \[\frac{\sqrt{n}(e^{-\overline{X}}-e^{-\lambda})}{\sqrt{\overline{X}e^{-2\overline{X}}}}\cong N(0,1).\] El IC \((1-\alpha)\%\) de Wald viene dado por: \(P_{\lambda}(\frac{\sqrt{n}\vert e^{-\overline{X}}-e^{-\lambda} \vert}{\sqrt{\overline{X}e^{-2\overline{X}}}}\leq qnorm(1-\alpha/2))\simeq 1-\alpha.\) De forma que despejando \(\lambda\) de la expresión anterior se llega al resultado:\[(e^{-\overline{X}}\pm qnorm(1-\alpha/2)\sqrt{\frac{\overline{X}e^{-2\overline{X}}}{n}})\]

g<-function(x){
  return(exp(-x))
}
g_prima<-function(x){
  return(-exp(-x))
}
emv_g<- g(emv)
emv_g

## [1] 0.2231302

emv.var_g<-emv.var*(g_prima(emv))^2
emv.var_g

##             [,1]
## [1,] 0.001493611

emv_g-qnorm(0.975)*sqrt(emv.var_g)

##           [,1]
## [1,] 0.1473829

emv_g+qnorm(0.975)*sqrt(emv.var_g)

##           [,1]
## [1,] 0.2988774

1. Calcular el IC basado en la RV para \(g(\lambda)=P_\lambda(X=0)\) con confianza aproximada del \(95\%.\)

Los IC basados en el estadístico de la RV son invariantes frente a transformaciones. Para calcular un IC para una función del parámetro bastará con aplicar dicha función a los extremos del intervalo RV obtenidos para el parámetro. En este caso el IC pedido es: \[(e^{-1.19},e^{-1.86})\simeq (0.16,0.30)\]

Ejercicio 2: Modelo de Exponencial.

Considere \(X_1,\dots,X_{30}\) m.a.s. de \(\mathcal{Exp}(\lambda),\) \(f(x,\lambda)=\lambda e^{-\lambda x},x>0\), \(\lambda>0.\) Para los datos observados se tiene \(\overline{X}=0.9\).

1. Obtener y representar la función de verosimilitud y log verosimilitud.
2. Calcular el EMV de forma analítica.
3. Obetener el EMV de forma numérica (aproximada) mediante la función {optim}.
4. Calcular IC de Wald para \(\lambda\) con confianza aproximada del \(95\%.\)
5. Calcular el estadístico test de Wald para contrastar \(H_0:\lambda=1 \quad \textit{vs}\quad H_1:\lambda\neq 1\). ¿Cuál es el p-valor del test?
6. Calcular el estadístico test de la RV para contrastar \(H_0:\lambda=1 \quad \textit{vs}\quad H_1:\lambda\neq 1\). ¿Cuál es el p-valor del test?
7. Calcular IC basado en la RV para \(\lambda\) con confianza aproximada del \(95\%.\)

Ejercicio 3:

Se seleccionan aleatoriamente 1000 muestras de suelo que se combinan en lotes de tamaño 10 muestras. Los lotes son examinados para evaluar la presencia de una toxina, es decir, las muestras individuales de suelo no se examinan. Concretamente, en 100 lotes, de 10 muestras de suelo cada uno, se detecto la toxina en 12 lotes.

1. Calcular el EMV para \(p=\text{probabilidad de que la toxina esté presente en una muestra de suelo individual}.\)
2. Detaerminar la distribución asintótica del EMV de \(p\) y calcular IC aproximados de Wald para \(p\) con confianza 0.95.
3. Calcular ICRV con confianza aproximada de 0.95 para \(p\).