Ejercicio 3.6

Consideremos \(N\) variables aleatorias binarias independientes \(Y_1, \dots, Y_N\) con \[ P(Y_i = 1) = \pi_i \quad \text{y} \quad P(Y_i = 0) = 1 - \pi_i \] La función de probabilidad de \(Y_i\) se puede escribir como \[ P(Y_i = y_i) = \pi_i^{y_i} (1 - \pi_i)^{1 - y_i} \] donde \(y_i = 0\) o \(1\).

(a) Demuestra que esta función de probabilidad pertenece a la familia exponencial de distribuciones.

Consideremos una única variable aleatoria \(Y\) cuya distribución de probabilidad depende de un único parámetro \(\theta\). La distribución pertenece a la familia exponencial si puede escribirse de la forma:\(f(y; \theta) = \exp\left[a(y)b(\theta) + c(\theta) + d(y)\right]\).

\[ f(y;\pi_i) = \pi_i^{y_i} (1 - \pi_i)^{1 - y_i} \\ = \exp\left\{ \ln \left( \pi_i^{y_i}(1-\pi_i)^{1-y_i} \right) \right\} \\ = \exp\left\{ y_i \ln \left(\pi_i \right) + (1 - y_i) \ln \left(1 - \pi_i \right) \right\} \\ = \exp\left\{ y_i \ln \left(\pi_i \right) + \ln \left(1 - \pi_i \right) - y_i \ln \left(1 - \pi_i \right) \right\} \\ = \exp\left\{ y_i \left[ \ln \left(\pi_i \right) - \ln \left(1 - \pi_i \right) \right] + \ln \left(1 - \pi_i \right) \right\} \\ = \exp\left\{ y_i \ln \left( \frac{\pi_i}{1 - \pi_i} \right) + \ln \left(1 - \pi_i \right) \right\} \]

Donde:

\[ a(y_i) = y_i\\ b(\pi_i) = \ln \left( \frac{\pi_i}{1 - \pi_i} \right)\\ c(\pi_i) = \ln(1 - \pi_i)\\ \]

no hay ningún término adicional que dependa solo de \(y_i\)

\[ d(y_i) = 0 \]

(b) De la demostración anterior podemos notar que el parámetro natural es:

\[ b(\pi_i) = \ln \left( \frac{\pi_i}{1 - \pi_i} \right) \]

Esta función, el logaritmo de los odds (razón de probabilidades” o simplemente “cociente de probabilidades”) \(\frac{\pi_i}{1 - \pi_i}\), se llama la función logit.

(c) Demuestra que \(E(Y_i) = \pi_i\)

El siguiente resultado pueden utilizarse ahora para distribuciones de la familia exponencial. \(E[a(Y)]=\frac{-c'(\theta)}{b'(\theta)}\)

Para la distribución binomial, sabemos que:

\[ E[a(Y_i)] =E[Y_i]= \frac{-c'(\pi_i)}{b'(\pi_i)} \]

Donde \(b(\pi_i) = \log\left(\frac{\pi_i}{1 - \pi_i}\right)\) y \(c(\pi_i) = \log(1 - \pi_i)\), se tiene que:

\[ -c'(\pi_i) = -\frac{d}{d\pi_i} \ln(1 - \pi_i) = \frac{1}{1 - \pi_i} \] Por otro lado: \[ b'(\pi_i) = \frac{1}{\frac{\pi_i}{1 - \pi_i}} \cdot \frac{d}{d\pi_i} \left( \frac{\pi_i}{1 - \pi_i} \right) \] Luego, usamos la regla del cociente para derivar:

\[ \frac{d}{d\pi_i} \left( \frac{\pi_i}{1 - \pi_i} \right) = \frac{(1 - \pi_i)(1) - \pi_i(-1)}{(1 - \pi_i)^2} = \frac{1}{(1 - \pi_i)^2} \]

Finalmente: \[ b'(\pi_i) = \frac{1}{\pi_i (1 - \pi_i)} \] Por lo tanto, la fórmula para la esperanza se convierte en:

\[ E[Y_i] = \frac{\frac{1}{(1 - \pi_i)}}{\frac{1}{\pi_i (1 - \pi_i)}} = \pi_i \]

(d) Si la función de enlace es

\[ g(\pi) = \log \left( \frac{\pi}{1 - \pi} \right) = \mathbf{x}^T \boldsymbol{\beta} \]

demuestra que esto es equivalente a modelar la probabilidad \(\pi\) como

\[ \pi = \frac{e^{\mathbf{x}^T \boldsymbol{\beta}}}{1 + e^{\mathbf{x}^T \boldsymbol{\beta}}} \]

Queremos expresar \(\pi\) en términos de \(\mathbf{x}^T \boldsymbol{\beta}\).

Despejamos \(\frac{\pi}{1 - \pi}\) de la ecuación:

\[ \log \left( \frac{\pi}{1 - \pi} \right) = \mathbf{x}^T \boldsymbol{\beta} \]

Tomamos la exponencial de ambos lados:

\[ \frac{\pi}{1 - \pi} = e^{\mathbf{x}^T \boldsymbol{\beta}} \] Multiplicamos ambos lados por \((1 - \pi)\):

\[ \pi = e^{\mathbf{x}^T \boldsymbol{\beta}} (1 - \pi) \]

Simplificamos:

\[ \pi = e^{\mathbf{x}^T \boldsymbol{\beta}} - e^{\mathbf{x}^T \boldsymbol{\beta}} \pi \]

\[ \pi + e^{\mathbf{x}^T \boldsymbol{\beta}} \pi = e^{\mathbf{x}^T \boldsymbol{\beta}} \]

\[ \pi \left( 1 + e^{\mathbf{x}^T \boldsymbol{\beta}} \right) = e^{\mathbf{x}^T \boldsymbol{\beta}} \]

\[ \pi = \frac{e^{\mathbf{x}^T \boldsymbol{\beta}}}{1 + e^{\mathbf{x}^T \boldsymbol{\beta}}} \]

(e) En el caso particular donde \(\mathbf{x}^T \boldsymbol{\beta} = \beta_1 + \beta_2 x\), esto da como resultado:

\[ \pi = \frac{e^{\beta_1 + \beta_2 x}}{1 + e^{\beta_1 + \beta_2 x}} \]

Para este caso tenemos que:

Entonces, la combinación lineal \(\mathbf{x}^T \boldsymbol{\beta}\) se simplifica a:

\[ \mathbf{x}^T \boldsymbol{\beta} = \beta_1 + \beta_2 x \] Resolviendo \(\log \left( \frac{\pi}{1 - \pi} \right) = \beta_1 + \beta_2 x\) obtenemos:

\[ \pi = \frac{e^{\beta_1 + \beta_2 x}}{1 + e^{\beta_1 + \beta_2 x}} \]

Esto es conocido como la función logística (Logit).

(f) Esboza el gráfico de \(\pi\) en función de \(x\) en este caso, tomando \(\beta_1,\beta_2\) como constantes. ¿Cómo interpretarías este gráfico si \(x\) es la dosis de un insecticida y \(\pi\) es la probabilidad de que un insecto muera?

Para esbozar el gráfico de \(\pi\) en función de \(x\), utilizamos la fórmula de la función logística:

\[ \pi = \frac{e^{\beta_1 + \beta_2 x}}{1 + e^{\beta_1 + \beta_2 x}} \]

# Cargar librerías necesarias
library(ggplot2)

# Definir las constantes
beta1 <- 0  # Intercepto
beta2_values <- c(-1, 0, 1)  # Diferentes valores de beta2

# Crear un rango de valores para x
x <- seq(-10, 10, length.out = 100)

# Crear un dataframe para almacenar los datos
data <- data.frame()

# Calcular pi para cada valor de beta2
for (beta2 in beta2_values) {
  pi <- exp(beta1 + beta2 * x) / (1 + exp(beta1 + beta2 * x))
  data <- rbind(data, data.frame(x = x, pi = pi, beta2 = factor(beta2)))
}

# Crear el gráfico
ggplot(data, aes(x = x, y = pi, color = beta2)) +
  geom_line(size = 1.2) +
  labs(title = "Probabilidad de muerte (\u03C0) en función de la dosis de insecticida (x)",
       x = "Dosis de insecticida (x)",
       y = "Probabilidad de muerte (\u03C0)",
       color = "\u03B2_2") +
  theme_minimal() +
  scale_color_manual(values = c("blue", "green", "red"), labels = c("-1", "0", "1"))

Interpretación del Gráfico

  • Cuando \(\beta_2 > 0\): El gráfico de \(\pi\) contra \(x\) será una curva sigmoidea ascendente. Esto indica que a medida que la dosis de insecticida aumenta, la probabilidad de que un insecto muera también aumenta.

  • Cuando \(\beta_2 < 0\): El gráfico será una curva sigmoidea descendente. Esto sugiere que a medida que la dosis de insecticida aumenta, la probabilidad de muerte del insecto disminuye. Aunque esto es menos común, podría ocurrir en situaciones particulares.

Este análisis es útil para comprender cómo varía la probabilidad de un evento (en este caso, la muerte de un insecto) en función de una variable predictora (la dosis de insecticida).

Ejercicio 3.7

¿Es la distribución de valor extremo (Gumbel), con función de densidad de probabilidad

\[ f(y; \theta) = \frac{1}{\phi} \exp \left\{ \frac{(y - \theta)}{\phi} - \exp \left[\frac{y - \theta}{\phi} \right] \right\} \]

(donde \(\phi > 0\) se considera un parámetro molesto) un miembro de la familia exponencial?.

\[ f(y; \theta) = \exp \left[ a(y) b(\theta) + c(\theta) + d(y) \right] \] Podemos reescribir como:

\[ f(y; \theta) = \exp \left\{ \ln \left[\ \frac{1}{\phi} \exp \left( \frac{(y - \theta)}{\phi} - \exp \left(\frac{y - \theta}{\phi} \right) \right) \right] \right\}\\ = \exp \left\{ \ln(1)-\ln(\phi) + \frac{(y - \theta)}{\phi} - \exp \left(\frac{y - \theta}{\phi} \right) \right\}\\ = \exp \left\{-\ln(\phi) + \frac{y}{\phi} - \frac{\theta}{\phi} - e^{\frac{y}{\phi}}e^{\frac{-\theta}{\phi}}\right\}\\ = \exp \left\{ - e^{\frac{y}{\phi}}e^{-\frac{\theta}{\phi}} - \frac{\theta}{\phi} -\ln(\phi) + \frac{y}{\phi} \right\} \]

Esta forma coincide con la forma general de la familia exponencial, donde:

Entonces, la distribución de valor extremo (Gumbel) es un miembro de la familia exponencial.

Ejercicio 3.8

Supongamos que \(Y_1, ..., Y_N\) son variables aleatorias independientes, cada una con la distribución de Pareto, y que

\[ E(Y_i) = (\beta_0 + \beta_1 x_i)^2. \]

¿Es este un modelo lineal generalizado? Da razones para tu respuesta.

Verifiquemos que el modelo lineal generalizado este definido en términos de:

1. Las distribuciones de todos los \(Y_i\) son de la misma forma (por ejemplo, todas Normales o todas Binomiales).

2. La distribución de cada \(Y_i\) tiene la forma canónica y depende de un único parámetro \(\theta_i\) (los \(\theta_i\) no tienen por qué ser todos iguales), por lo que

Nota: Si \(a(y) = y\), se dice que la distribución está en forma canónica (es decir, estándar) y \(b(\theta)\) se denomina a veces parámetro natural de la distribución.

\[ f(y_i; \theta_i) = \exp \left[ y_i b(\theta_i) + c(\theta_i) + d(y_i) \right]. \] Por el enunciado, sabemos que contamos con un conjunto de variables aleatorias continuas independientes, \(Y_i\sim Pareto\left ( \theta_i , 1 \right )\) con \((i=1,...N)\) Y función de densidad de probabilidad es:

\[ f(y_i; \theta) = \theta_i y_i^{-\theta_i - 1}= \theta_i y_i^{-(\theta_i + 1)} \]

Veamos si las \(Y_i\) pertenecen a la famila exponencial:

\[ f(y_i; \theta_i) = \exp \left\{\ln \left [ \theta_i y_i^{-(\theta_i + 1)} \right ]\right\}\\ = \exp \left\{\ln (\theta_i) - (\theta_i + 1)\ln(y_i)\right\}\\ = \exp \left\{\ln (\theta_i) - \theta_i\ln(y_i) + \ln(y_i)\right\}\\ = \exp \left\{\theta_i \cdot \left [-\ln(y_i)\right ] + \ln (\theta_i) + \ln(y_i)\right\} \] Hemos demostrado que la distribución pareto sí pertenece a la familia exponencial de distribuciones donde:

Con los resultados anteriores podemos notar que no se cumple que la distribución está en forma canónica, ya que \(a(y_i) \neq y_i\). Por lo que podemos decir que el modelo lineal generalizado no está definido para esta distribución, por lo que no se trata de un modelo lineal generalizado.

Sabemos que \(\mu_i = E(Y_i)\) es la media de la distribución.

En este caso, tenemos:

\[ \mu_i=E(Y_i) = (\beta_0 + \beta_1 x_i)^2 \] En este caso, sí para la distribución pareto el GLM estuviera definido, la función de enlace debe ser la raíz cuadrada:

\[ g(\mu_i) =\sqrt{ E(Y_i)} = \beta_0 + \beta_1 x_i \]

Ejercicio 3.9

Sean \(Y_1, \dots, Y_N\) variables aleatorias independientes con:

\[ E(Y_i) = \mu_i = \beta_0 + \log (\beta_1 + \beta_2 x_i) \quad ; \quad Y_i \sim N(\mu, \sigma^2) \]

para todos \(i = 1, \dots, N\).

Para determinar si este es un modelo lineal generalizado (GLM), debemos revisar los componentes básicos de un GLM:

  1. Distribución de la variable dependiente: Un GLM debe tener una distribución perteneciente a la familia exponencial. En este caso, la variable dependiente \(Y_i\) sigue una distribución Normal, que pertenece a la familia exponencial, lo que es consistente con un GLM.(se asume que \(\sigma^{2}\) es conocido para reescribir \(f(y; \theta) = \exp \left[ a(y) b(\theta) + c(\theta) + d(y) \right]\))

\[ f(y; \mu) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( -\frac{(y - \mu)^2}{2\sigma^2} \right)\\ = \exp \left\{\ln \left [ \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( -\frac{(y - \mu)^2}{2\sigma^2} \right)\right ] \right\}\\ = \exp \left\{\ln(1)-\ln(\sqrt{2\pi \sigma^2}) -\frac{(y - \mu)^2}{2\sigma^2} \right\}\\ = \exp \left\{ -\frac{1}{2} \ln(2\pi \sigma^2) -\frac{y^{2}}{2\sigma^{2}} + \frac{2y\mu}{2\sigma^{2}} - \frac{\mu^{2}}{2\sigma^{2}} \right\}\\ = \exp \left\{ y\frac{\mu}{\sigma^{2}} - \frac{\mu^{2}}{2\sigma^{2}}-\frac{1}{2} \ln(2\pi \sigma^2) -\frac{y^{2}}{2\sigma^{2}} \right\} \]

Podemos identificar que:

Con los resultados anteriores podemos observar que se cumple que, la distribución está en forma canónica, ya que \(a(y)= y\).

  1. Función de enlace: En un GLM, existe una función de enlace \(g(\mu_i)\) que conecta el valor esperado de la variable dependiente \(\mu_i = E(Y_i)\) con una combinación lineal de las covariables.

En este caso, el valor esperado es \(\mu_i = \beta_0 + \log(\beta_1 + \beta_2 x_i)\) Si se toma la función exponencial como función de enlace, entonces…

\[ g(\mu_i)= \exp[E(Y_i)] = \exp[\beta_0 +log(\beta_1 +\beta_2x_i)]\\ =e^{\beta_0} e^{log(\beta_1 +\beta_2 x_i)} \\ = e^{\beta_0} (\beta_1 +\beta_2 x_i)\\ = e^{\beta_0} \beta_1+ e^{\beta_0}\beta_2 x_i\\ \] En este contexto, el término \(e^{\beta_0}\) actúa como un factor de escala para los parámetros \(\beta_1\) y \(\beta_2\), ajustando su impacto en la función de enlace. Este término debe ser una constante conocida.

Donde podemos definir nuevos parámetros:

\[ \beta_1^* = e^{\beta_0} \beta_1\\ \beta_2^* = e^{\beta_0} \beta_2 \]

Entonces, la forma simplificada es:

\[ g(\mu_i) = \beta_1^* + \beta_2^* x_i \] La función de enlace resultante es una combinación lineal de las covariables \(x_i\). Esta forma coincide con la estructura de los modelos lineales generalizados (GLM), donde el valor esperado de la variable dependiente está relacionado linealmente con las covariables a través de la función de enlace.

Por lo tanto, el modelo dado es un modelo lineal generalizado.

Ejercicio 3.10

Para la distribución Pareto, encuentre la estadística score \(U\) y la información \(\mathfrak{J} =I= \text{var}(U)\). Verifique que \(E(U) = 0\).

La estadística score \(U\) A partir de la definición de familia exponencial \(f(y; \theta) = \exp \left[ a(y) b(\theta) + c(\theta) + d(y) \right]\)), la función de log-verosimilitud para una distribución de la familia exponencial es:

\[ l(\theta; y) = a(y)b(\theta) + c(\theta) + d(y) \] La derivada de \(l(\theta; y)\) con respecto a \(\theta\) es:

\[ U(\theta; y) = \frac{d l(\theta; y)}{d \theta} = a(y)b'(\theta) + c'(\theta) \]

La función \(U\) se denomina estadístico de puntuación (estadística score) y, como depende de \(y\), puede considerarse una variable aleatoria, es decir:

\[ U = a(Y)b'(\theta) + c'(\theta) \]

En el ejercicio 3.8 Hemos demostrado que la distribución pareto sí pertenece a la familia exponencial de distribuciones donde:

Las derivadas correpondientes son \(b'(\theta)=-1\) y \(c'(\theta)=1/\theta\), remplazando obtenemos \(U\).

\[ U = -\ln(y)+\frac{1}{\theta} \]

Estos resultado pueden utilizarse para distribuciones de la familia exponencial.

Para el cálculo de \(\text{Var}[a(Y)]\) la primera y segundas derivadas correspondientes \(c''(\theta)=- \frac{1}{\theta^2}\), \(b''(\theta) = 0\), Remplazando obtenemos:

\[ \text{Var}[a(Y)] = \frac{b''(\theta) c'(\theta) - c''(\theta) b'(\theta)}{[b'(\theta)]^3}\\ = \frac{(0) c'(\theta) - (-\frac{1}{\theta^2}) (-1))}{[-1]^3}\\ = \frac{(-\frac{1}{\theta^2})}{[-1]}\\ = \frac{1}{\theta^2} =\theta^{-2} \]

El valor esperado de \(U\) es:

La Información \(I = \text{var}(U)\)

Por la fórmula de la varianza de una transformación lineal de variables aleatorias \(V(\alpha X+ \lambda)= \alpha^{2} V(X)\).

  • \(I = \text{var}(U) = \text{var}[a(Y)b'(\theta) + c'(\theta)]=\left[b'(\theta)\right]^2 \text{var}[a(Y)]\). Reemplzando los valores obtenidos anteriormente \(\mathfrak{J} =I= [-1]^{2}\theta^{-2}\)

Otra forma alternativa para calcular \(I\) sin tener que calcular la \(Var[a(Y)]\) es usar \(\text{var}(U) = \frac{b''(\theta) c'(\theta)}{b'(\theta)} - c''(\theta)=\frac{(0) c'(\theta)}{b'(\theta)} -(- \frac{1}{\theta^{2}})= \theta^{-2}\).

Ejercicio 4.1

(Cap 4, 4.1,4.3 - Cap 3, 3.6:3.10)

Los datos de la Tabla 4.5 (Número de casos de SIDA en Australia por trimestres sucesivos de 1984 a 1988) muestran el número de casos de SIDA en Australia por fecha de diagnóstico durante periodos sucesivos de 3 meses de 1984 a 1988. (Datos del National Centre for HIV Epidemiology and Clinical Research, 1994).

En esta primera fase de la epidemia, el número de casos parecía aumentar exponencialmente.

a. Represente gráficamente el número de casos \(y_{i}\) frente al periodo de tiempo i (i =1,…,20).

# Crear un vector con los números de casos
cases <- c(1, 6, 16, 23, 27, 39, 31, 30, 43, 51, 63, 70, 88, 97, 91, 104, 110, 113, 149, 159)

# Crear un vector con los periodos de tiempo (i)
time_period <- 1:20

plot(time_period, cases)

# Graficar los casos contra el tiempo
plot(time_period, cases, type="o", col="blue", xlab="Time Period (i)", ylab="Number of Cases", main="Number of AIDS Cases in Australia (1984-1988)")

# Agregar una línea que conecte los puntos
lines(time_period, cases, col="blue")

  1. Un modelo posible es la distribución de Poisson con el parámetro \(\lambda_i= i^\theta\), o equivalentemente \(\log \lambda_i = \theta \log i\). Grafique \(\log y_i\) contra \(\log i\) para examinar este modelo.
# Calcular los logaritmos de los casos y de los periodos de tiempo

log_time_period <- log(time_period)

# Graficar log(y_i) contra log(i)
plot(log_time_period, cases, type="o", col="red", xlab="log(i)", ylab="log(y_i)", main="Log Plot of AIDS Cases in Australia (1984-1988)")

plot(log_time_period, cases)

c. Ajusta un modelo lineal generalizado a estos datos utilizando la distribución de Poisson, la función de enlace logarítmica y la ecuación \(g(\lambda_i) = \log \lambda_i = \beta_1 + \beta_2 x_i\), donde \(x_i = \log i\). En primer lugar, haz esto desde los primeros principios, desarrollando las expresiones para la matriz de pesos \(W\) y otros términos necesarios para la ecuación iterativa \[X^T W X \mathbf{b}^{(m)} = X^T W \mathbf{z}\] y utilizando software que pueda realizar operaciones matriciales para llevar a cabo los cálculos.

Datos4.5 <- cbind(log_time_period, cases)

Datos4.5 <- as.data.frame(Datos4.5)
X<- matrix(c(rep(1,20), log(1:20)), ncol = 2)

W<-function(X,B){
  n<-dim(X)[1]
  I<- diag(c(exp(X%*%B)))
  return(I)
}

IF<- function(X,B){
  return(matrix(c(round(sum(X[,1]^2*exp(X%*%B)),3),
                  round(sum(X[,1]*X[,2]*exp(X%*%B)),3),
                  round(sum(X[,1]*X[,2]*exp(X%*%B)),3),
                  round(sum(X[,2]^2*exp(X%*%B)),3)), ncol=2))
}

z <- function(X,B,Y){
  return(X%*%B +(Y-exp(X%*%B))/exp(X%*%B))
}

B0<-c(0.2, 2)
B<-c(1, 4)

while(all(abs(B-B0) > c(0.002,0.002))== TRUE){
  B1 <- solve(t(X)%*%W(X,B0)%*%X)%*%t(X)%*%W(X,B0)%*%z(X,B0,Y=Datos4.5$cases)
  B <- B0
  B0 <- B1
}

B1
         [,1]
[1,] 0.995998
[2,] 1.326610

d. Ajusta el modelo descrito en (c) utilizando software estadístico que pueda realizar regresión de Poisson. Compara los resultados con los obtenidos en (c).

Datos4.5 <- cbind(log_time_period, cases)

Datos4.5 <- as.data.frame(Datos4.5)


fit <- glm(formula =cases~log_time_period, family=poisson(link="log"), data =Datos4.5)

summary(fit)  

Call:
glm(formula = cases ~ log_time_period, family = poisson(link = "log"), 
    data = Datos4.5)

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.99600    0.16971   5.869 4.39e-09 ***
log_time_period  1.32661    0.06463  20.525  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 677.264  on 19  degrees of freedom
Residual deviance:  21.755  on 18  degrees of freedom
AIC: 138.05

Number of Fisher Scoring iterations: 4

Graficamente

fun <- function(x){
  exp(coef(fit)[1]+coef(fit)[2]*x)
}

y1<-fun(log_time_period)
df <- data.frame(x=log_time_period,y=cases,y1)
library(ggplot2)
ggplot(df, aes(x)) +                    # basic graphical object
  geom_point(aes(y=y), colour="red") +  # first layer
  geom_line(aes(y=y1), colour="green")+geom_function(fun = fun, colour = "blue")+
  labs(x = "log10(initial white blood cell count)")+
  labs(x = "log(i)")+
  labs(title = "Log Plot of AIDS Cases in Australia (1984-1988)", subtitle = "MLG de La Paz")

Ejercicio 4.3

Sea \(Y_1, ..., Y_N\) una muestra aleatoria de la distribución Normal \(Y_i \sim N(\log \beta, \sigma^2)\), donde \(\sigma^2\) es conocida. Encontrar el estimador de máxima verosimilitud de \(\beta\). También verificar las ecuaciones \(U_j = \sum_{i=1}^{N} \left( \frac{(y_i - \mu_i)}{\text{var}(Y_i)} x_{ij} \frac{\partial \mu_i}{\partial \eta_i} \right)\) y \(\mathbf{X}^T \mathbf{W} \mathbf{X} \boldsymbol{b}^{(m)} = \mathbf{X}^T \mathbf{W} \mathbf{z}\) en este caso.

Estimación de máxima verosimilitud para \(\beta\)

\[ f(y_{1},..,y_{N};\beta)= \prod_{i=1}^{N} \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( -\frac{(y_i -log(\beta))^2}{2\sigma^2} \right)\\ =\left(\frac{1}{\sqrt{2\pi \sigma^2}}\right)^{N} \exp \left ( \frac{-1}{2\sigma^{2}}\sum_{i=1}^{N}(y_i -log(\beta))^2\right )\\ = \left (\sqrt{2\pi \sigma^2} \right )^{-N} \exp \left ( \frac{-1}{2\sigma^{2}}\sum_{i=1}^{N}(y_i -log(\beta))^2\right ) \]

El logarítmo de la función de verosimilitud es:

\[ l= -N \ln (\sigma\sqrt{2\pi})- \frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(y_i -log(\beta))^2 \]

Derivando e igualando a cero para determinar puntos críticos, Para encontrar el \(\beta\) que máximiza la función de log-verosimilitud.

\[ \frac{ \partial l }{\partial \beta} = - \frac{2}{2\sigma^{2}} \left (- \frac{1}{\beta} \right ) \sum_{i=1}^{N}(y_i -\log(\beta))\\ = \frac{1}{\sigma^{2}} \frac{1}{\beta} \sum_{i=1}^{N}(y_i -\log(\beta)) \]

Entonces: \[ \frac{1}{\sigma^{2}} \frac{1}{\beta} \sum_{i=1}^{N}(y_i -\log(\beta)) = 0\\ \sum_{i=1}^{N}(y_i) = N\log(\beta)\\ \log(\beta)= \bar{Y}\\ \beta= e^{\bar{Y}} \]

El estimador de máxima verosimilitud para nuestro parametro \(\hat{\beta}=\frac{1}{N}\sum_{i=1}^{N}Y_i\).

La estadística Score

\(U_j = \sum_{i=1}^{N} \left( \frac{(y_i - \mu_i)}{\text{var}(Y_i)} x_{ij} \frac{\partial \mu_i}{\partial \eta_i} \right)\).

Sabemos que: \(E(Y_i)=\mu_i=\log(\beta)\), para este caso la función de enlace es: \(g(\mu_i)=e^{E(Y_i)}= \beta=\eta_i\). Por lo tanto tenemos que \(U\) es:

\[ U = \sum_{i=1}^{N} \left( \frac{(y_i - \log(\beta))}{\sigma^2} \frac{1}{\beta} \right) \]

Donde \(x_{ij}=1\), \(\frac{\partial \mu_i}{\partial \eta_i}=\frac{1}{\beta}\)

Ejercicio 5.1

(a) Encuentra la estadística de Wald \((\widehat{\pi} - \pi)^T \mathfrak{J} (\widehat{\pi} - \pi)\), donde \(\widehat{\pi}\)es el estimador de máxima verosimilitud de \(\pi\) y \(\mathfrak{J}\) es la información.

Iniciamos calculando el estiamdor de máxima verosimilitud para \(\pi\) La función de probabilidad es: \(f(y;\pi) =\binom{n}{y} \pi^{y} (1 - \pi)^{n - y}\)

La función de log-verosimilitud que corresponde al logaritmo de la función de densidad conjunta, pero en este caso tenemos una unica variable aleatoria.

\[ l(\pi;y) = \ln \prod_{i=1}^{1} f(y;\pi)\\ = \ln [\binom{n}{y} \pi^{y} (1 - \pi)^{n - y}]\\ = \ln[\binom{n}{y}]+ y\ln(\pi)+(n-y)\ln(1-\pi)\\ = \ln[\binom{n}{y}]+ y\ln(\pi)+n\ln(1-\pi)-y\ln(1-\pi) \]

Derivamos la función de log-verosimilitud \(U(\pi;y)=\frac{d l(\pi; y)}{d \pi}= \frac{Y}{\pi}-\frac{n}{(1-\pi)}+\frac{y}{(1-\pi)}\). Igualando a cero y depejando \(\pi\) para deteminar el estiamdor de máxima verosimilitud.

\[ U=\frac{\partial l(\pi;y)}{\partial \pi}=\frac{Y}{\pi}-\frac{n}{(1-\pi)}+\frac{y}{(1-\pi)}=0\\ \frac{Y}{\pi}+\frac{Y-n}{(1-\pi)}=0\\ \frac{Y}{\pi}=-\frac{Y-n}{(1-\pi)}\\ {Y}(1-\pi)=-y\pi+n\pi\\ Y-{Y}\pi=-Y\pi+n\pi\\ \text{El estimador de máxima verosimilitud es: }\widehat{\pi}=\frac{Y}{n} \\ \text{ De modo que las estimación de máxima verosimilitud }\widehat{\pi}=\frac{y}{n} \]

Ahora calculamos \(Var(U)=\mathfrak{J}\) que esta dado por:\(U = \frac{dl}{d\pi} = \frac{Y}{\pi} - \frac{n - Y}{1 - \pi} = \frac{Y - n\pi}{\pi(1 - \pi)}\)

\[ Var(U)=I=\frac{Var(Y)}{(\pi(1 - \pi))^{2}}\\ = \frac{n\pi(1 - \pi)}{(\pi(1 - \pi))^{2}}= \frac{n}{\pi(1 - \pi)} \]

Entonces la estadística de Wald esta dada como:

\[ (\widehat{\pi} - \pi)^T \mathfrak{J} (\widehat{\pi} - \pi) = \left(\frac{y}{n} - \pi\right) \frac{n}{\pi(1 - \pi)} \left(\frac{y}{n} - \pi\right) \\ = \left(\frac{y}{n} - \pi\right)^2 \frac{n}{\pi(1 - \pi)} = \frac{(y - n\pi)^2}{n \pi(1 - \pi)} \]

(b) Verifica que la estadística de Wald es la misma que la estadística \(U^T \mathfrak{J}^{-1}U\) en este caso.

Debemos mostrar que \[(\widehat{\pi} - \pi)^T \mathfrak{J} (\widehat{\pi} - \pi)=U^T \mathfrak{J}^{-1}U\] sabemos que \[(\widehat{\pi} - \pi)^T \mathfrak{J} (\widehat{\pi} - \pi)=\frac{(y - n\pi)^2}{n \pi(1 - \pi)}\], \[\mathfrak{J}=\frac{n}{\pi(1 - \pi)}\], \[U=\frac{Y - n\pi}{\pi(1 - \pi)}\]. Reemplazando obtenemos:

\[ U^T \mathfrak{J}^{-1}U= \left ( \frac{Y - n\pi}{\pi(1 - \pi)} \right ) \left ( \frac{n}{\pi(1 - \pi)} \right )^{-1} \left ( \frac{Y - n\pi}{\pi(1 - \pi)} \right ) \\ = \left ( \frac{Y - n\pi}{\pi(1 - \pi)} \right )^{2} \frac{\pi(1 - \pi)}{n}\\ = \frac{(y - n\pi)^2}{n \pi(1 - \pi)} \]

(c) Encuentra la desviación \[ 2[l(\widehat{\pi}; y) - l(\pi; y)] \] Tenemos que: \[ l(\widehat{\pi};y) = \ln[\binom{n}{y}]+ y\ln(\widehat{\pi})+(n-y)\ln(1-\widehat{\pi})\\ l(\pi;y)=\ln[\binom{n}{y}]+ y\ln(\pi)+(n-y)\ln(1-\pi)\\ D= 2\left [\ y\ln(\widehat{\pi})+(n-y)\ln(1-\widehat{\pi}) - y\ln(\pi)- (n-y)\ln(1-\pi) \right] \\ \text{Agrupamos los términos que contienen \( y \) y los que contienen \( n - y \)}\\ D = 2 \left[ y \left( \ln(\widehat{\pi}) - \ln(\pi) \right) + (n - y) \left( \ln(1 - \widehat{\pi}) - \ln(1 - \pi) \right) \right]\\ \text{Usamos la propiedad de los logaritmos \( \ln(a) - \ln(b) = \ln\left(\frac{a}{b}\right) \) para simplificar:}\\ D = 2 \left[ y \ln\left(\frac{\widehat{\pi}}{\pi}\right) + (n - y) \ln\left(\frac{1 - \widehat{\pi}}{1 - \pi}\right) \right] \]

Donde:\(\widehat{\pi}=\frac{y}{n}\)

(d) Para muestras grandes, tanto la estadística de Wald como la desviación tienen aproximadamente la distribución \(\chi^2(1)\). Para \(n = 10\) y \(y = 3\), usa ambas estadísticas para evaluar la adecuación de los modelos: 1. \(\pi = 0.1\), 2. \(\pi = 0.3\), 3. \(\pi = 0.5\).

¿Llevan ambas estadísticas a las mismas conclusiones?

La distribución \(\chi^2(1)\) tiene un percentil 95 de 3.84, el cual puede ser utilizado como el valor crítico:

# Parámetros
n <- 10
y <- 3
pi_values <- c(0.1, 0.3, 0.5)

# Inicializar vectores para almacenar resultados
wald_stats <- numeric(length(pi_values))
deviance_stats <- numeric(length(pi_values))

# Calcular estadísticas para cada valor de pi
for (i in 1:length(pi_values)) {
  pi <- pi_values[i]
  
  # Estadística de Wald
  wald_stats[i] <- (y - n * pi)^2 / (n * pi * (1 - pi))
  
  # Estadística de desviación
  pi_E <- y / n
  deviance_stats[i] <- 2 * (y * log(pi_E / pi) + (n - y) * log((1 - pi_E) / (1 - pi)))
}

# Crear tabla de resultados
tabla <- data.frame(
  Pi = pi_values,
  Wald_Statistic = wald_stats,
  Deviance_Statistic = deviance_stats
)

# Mostrar resultados
print(tabla)
   Pi Wald_Statistic Deviance_Statistic
1 0.1       4.444444           3.073272
2 0.3       0.000000           0.000000
3 0.5       1.600000           1.645658

Los resultados indican que la estadística de Wald/score sugiere rechazar \(\pi = 0.1\) (ya que 4.44 es mayor que el valor crítico de 3.84), mientras que la estadística de log-verosimilitud no sugiere rechazarlo (ya que 3.07 es menor que el valor crítico de 3.84).

Para \(\pi = 0.3\), ambas estadísticas son cero, lo que no sugiere rechazar \(\pi = 0.3\), ya que cero es mucho menor que el valor crítico de 3.84.

Tanto la estadística de Wald/score como la de log-verosimilitud para \(\pi = 0.5\) son menores que el valor crítico de 3.84, por lo que ninguna de las dos sugiere rechazar \(\pi = 0.5\).

Ejercicio 5.2

\[f(y_i; \theta_i) = \theta_i \exp(-y_i \theta_i)\]

Deriva la desviación comparando el modelo máximo con diferentes valores de \(\theta_i\) para cada \(Y_i\) y el modelo con \(\theta_i = \theta\) para todos los \(i\).

\(D = 2 \left[ l({\theta}_{\text{máx}}; y_{i}) - l(\theta; y) \right]\) Para calcular en desviación debemos calcular la log_verosimilitud para encontrar el modelo máximo o saturado.

\[ l({\theta}_{\text{máx}};y_{i}) = \ln \prod_{i=1}^{N} f(y_i; \theta_i)\\ \text{Usando la propiedad del logaritmo del producto, se simplifica a:} \\ = \sum_{i=1}^{N} \ln f(y_i; \theta_i)\\ = \sum_{i=1}^{N} \ln \left [ \theta_i \exp(-y_i \theta_i) \right ]\\ = \sum_{i=1}^{N} \left [ \ln(\theta_{i}) - y_i \theta_i\right ] \]

Derivando e igualando a cero, para encontarar un punto crítico tal que maximize la función Log_verosimilitud. \[ \frac{\partial l({\theta}_{\text{máx}};y_{i})}{\partial \theta_{i}}= \frac{1}{\theta_{i}} - y_{i} = 0\\ \frac{1}{\theta_{i}}=y_{i} \\ \theta_{i}= \frac{1}{y_{i}}\\ \text{La estimación de maxima verosimilitud es}\\ \widehat{\theta_{i}}= \frac{1}{y_{i}} \]

De modo que el estimador de maxima verosimilitud correponde a:\(\widehat{\theta_{i}}= \frac{1}{Y_{i}}\)

Por otro lado: para \(\theta_{i}= \theta\) \[ l(\theta; y) = \sum_{i=1}^{N} \left [ \ln(\theta) - y_{i} \theta \right ]\\ = N\ln(\theta)- \theta \sum_{i=1}^{N} y_{i} \]

Derivando e igualando a cero, para encontarar un punto crítico tal que maximize la función Log_verosimilitud.

\[ \frac{\partial l(\theta; y)}{\partial \theta} = \frac{N}{\theta} - \sum_{i=1}^{N} y_{i} = 0 \\ \frac{N}{\theta} = \sum_{i=1}^{N} y_{i} \\ \text{La estimación de máxima verosimilitud es} \\ \widehat{\theta} = \frac{N}{\sum_{i=1}^{N} y_{i}} = \frac{1}{\frac{\sum_{i=1}^{N} y_{i}}{N}} \\ \widehat{\theta} = \frac{1}{\bar{y}} \] De modo que el estimador de maxima verosimilitud correponde a:\(\widehat{\theta} = \frac{1}{\bar{Y}}\)

Entonces el desvio esta dado por:

\[ D=2\left[ l({\theta}_{\text{máx}}; y_{i}) - l(\theta; y) \right]\\ = 2\left[ \sum_{i=1}^{N} \left [ \ln(\widehat{\theta_{i}}) - y_i \widehat{\theta_{i}}\right ]- \left [ N\ln(\theta)- \theta \sum_{i=1}^{N} y_{i} \right] \right]\\ \text{Ya hemos calulado } \widehat{\theta_{i}}, \text{Sabemos que } \widehat{y_{i}}=E[y_{i}]=\frac{1}{\theta}, \text{entonces }\theta=\frac{1}{\widehat{y_{i}}} \\ = 2 \sum_{i=1}^{N} \left[ \ln\left(\frac{1}{y_{i}}\right) - \frac{y_{i}}{y_{i}} - \ln\left(\frac{1}{\widehat{y_{i}}}\right) + \frac{y_{i}}{\widehat{y_{i}}} \right]\\ = 2 \sum_{i=1}^{N} \left[ -\ln(y_{i}) - 1 - (-\ln(\widehat{y_{i}})) + \frac{y_{i}}{\widehat{y_{i}}} \right]\\ = 2 \sum_{i=1}^{N} \left[ \ln(\widehat{y_{i}}) - \ln(y_{i}) - 1 + \frac{y_{i}}{\widehat{y_{i}}} \right]\\ = 2 \sum_{i=1}^{N} \left[ -\ln\left(\frac{y_{i}}{\widehat{y_{i}}}\right) - 1 + \frac{y_{i}}{\widehat{y_{i}}} \right] \]

5.3

Supongamos que \(Y_1, \ldots, Y_N\) son variables aleatorias independientes e idénticamente distribuidas con la distribución de Pareto con parámetro \(\theta\).

(a)Encuentra el estimador de máxima verosimilitud \(\hat{\theta}\) de \(\theta\).

La función de densidad para la distribución pareto con parametros \((\theta,1)\):

\[f(y_i; \theta) = \theta y_i^{-\theta - 1}\]

\[f(\theta; y_1, \ldots, y_N) = \prod_{i=1}^{N} f(y_i; \theta) = \prod_{i=1}^{N} \theta y_i^{-\theta - 1}\]

La log_verosimilitud es:

$$ l(; y_1, , y_N) = {i=1}^{N} (y_i^{-- 1})\ l(; y_1, , y_N, 1) = {i=1}^{N} \ \ = _{i=1}^{N} ( - (y_i))\

\ = _{i=1}^{N} (y_i) \ = $$

(b)Encuentra la estadística de Wald para hacer inferencias sobre \(\theta\) (Sugerencia: Usa los resultados del Ejercicio 3.10).

\[ (\hat{\theta} - \theta)^T \mathfrak{J}(\theta) (\hat{\theta} - \theta) \sim \chi^2(1) \]

Se puede reescribir como: \[ (\hat{\theta} - \theta)^T \mathfrak{J}(\theta) (\hat{\theta} - \theta)\\ = (\hat{\theta} - \theta)^2 \mathfrak{J}(\theta) \]

En el ejercicio 3.10 y por el ejercicio 3.8, se obtuvo que \((\mathfrak{J}(\theta) = \text{Var}(a(Y_i))\) \(a(Y_i) = \sum_{i=1}^{N} \ln(Y_i)\), Entonces:

\[ \mathfrak{J}(\theta) = \text{Var}(a(Y_i)) = \text{Var} \left( \sum_{i=1}^{N} \ln(Y_i) \right)\\ \text{Por independencia}\\ \mathfrak{J}(\theta) = \text{Var}(a(Y_i)) = \sum_{i=1}^{N} \text{Var}(\ln(Y_i))\\ \text{Como se demostró en el ejercico 3.10}\\ \text{Var}(\ln(Y_i)) = \frac{1}{\theta^2} \]

Por lo tanto

\[ \mathfrak{J}(\theta) = \text{Var}(a(Y_i)) = \sum_{i=1}^{N} \frac{1}{\theta^2}\\ = \frac{N}{\theta^2} \] Reemplzando las expresiones en el estadistico de wald \[ (\theta-\hat{\theta})^2 \frac{1}{N \theta^2}\\ = \left( \theta - \frac{N}{\sum_{i=1}^{N} \ln(y_i)} \right)^2 \frac{1}{N \theta^2} \]

(c)Usa la estadística de Wald para obtener una expresión para un intervalo de confianza aproximado del 95% para \(\theta\).

(d)Las variables aleatorias \(Y\) con la distribución de Pareto con el parámetro \(\theta\) se pueden generar a partir de números aleatorios \(U\) que están uniformemente distribuidos entre 0 y 1 usando la relación \[ Y = \left(\frac{1}{U}\right)^{1/\theta} \] (Evans et al., 2000). Usa esta relación para generar una muestra de 100 valores de \(Y\) con \(\theta = 2\). A partir de estos datos, calcula una estimación \(\hat{\theta}\). Repite este proceso 20 veces y también calcula intervalos de confianza del 95% para \(\theta\). Compara el promedio de las estimaciones \(\hat{\theta}\) con \(\theta = 2\). ¿Cuántos de los intervalos de confianza contienen \(\theta\)?