En el fichero Datos-Encuesta-USA.RData aparecen datos del General Social Survey hecho en USA a finales de los 90, hemos seleccionado sólo 5 de las variables presentes en la encuesta, concretamente:

Lee los datos con ayuda de R.

load("DatosEncuestaUSA.RData")

El objetivo es analizar el número de hijos que tienen las mujeres que tenían una edad de 40 años en el instante de participación en la encuesta, interesa analizar las mujeres que participaron en la encuesta en la década de 1990.

Con el siguiente código R puedes seleccionar la variable de interés, hijos, de mujeres con 40 años y que particiaron en los años 90 en la encuesta (la encuesta no incluye años más allá de 1999).

hijos<-Y$CHILDS[Y$FEMALE==1 & Y$YEAR>=1990 & Y$AGE==40 & !is.na(Y$DEG)]
hijos<-hijos[!is.na(hijos)]

En el objeto hijos aparecen el número de hijos que tienen las mujeres seleccionadas. En concreto se trata de una muestra de tamaño \(n = 155\).

Es muy usual en la práctica de la estadística usar el modelo Poisson para representar una variable aleatoria que es un conteo, veremos más adelante cómo usar técnicas gráficas y procedimientos para analizar suposiciones como ésta (asumir un determinado modelo paramétrico para uos datos). Por el momento, supongamos que \[N_i = \mbox{ Número de hijos que tiene cada mujer de esta población}\]

esta variable es aleatoria, y representamos su aleatoriedad mediante una \(Po(\lambda)\), con \(\lambda\) el número esperado de hijos que tienen las mujeres de esta poblacion.

Es decir la masa de probabilidad, suponiendo \(\lambda\) conocido es:

\[P(X = x) = \frac{e^{-\lambda} \lambda^{x}}{x!}\ ; x = 0, 1, 2,..., \infty\]

El objetivo de este ejercicio es estimar puntualmente \(\lambda\).

Estimador máximo verosímil

Escribe la función de verosimilitud para \(\lambda\). Recuerda que es:

\[L(\lambda | n_1,\ldots,n_{155}) = f(n_1, \ldots, n_{155} | \lambda)\]

Escribe esta función para el caso del modelo Poisson. Encuenta el valor de \(\lambda\) que lo maximiza. Puedes usar técnicas matemáticas para encontrar el máximo, o calcularlo con ayuda de R (dibujando la función, donde el argumento es \(\lambda\), y usando la función optimize de R para buscar el valor que maximiza esta función).

Para buscar el máximo usando R tienes que programar primero la función verosimilitud en R (yo la he llamado vero.pois), y después usar los siguientes argumentos:

optimize(f = , interval = , …, lower = min(interval), upper = max(interval), maximum = FALSE, tol = .Machine$double.eps^0.25)

Tienes que especificar que optimize calcule un máximo, y que busque en un intervalo (por ejemplo prueba con el intervalo 0, 10), además hay que pasarle el resto de argumentos de la función vero.pois.

Desarrollo del ejercicio

\[ L(\lambda | x_i) = \prod_{i = 1}^{155} \frac{e^{-\lambda} · \lambda^{x_i}}{x_i!} = \frac{e^{-n·\lambda}·\lambda^{\sum x_i}}{\prod x_i!} \]

vero.poiss <- function(lambda) {
  n <- length(hijos)
  return((exp(-n*lambda)*lambda^(sum(hijos))) / (prod(factorial(hijos))))
}

interval <- 0:10

optimize(f = vero.poiss, 
         interval = interval, 
         lower = min(interval), 
         upper = max(interval), 
         maximum = TRUE
         )
## $maximum
## [1] 1.825806
## 
## $objective
## [1] 1.941674e-112

\(\hat{\lambda} = 1.825806\)

Vamos a estimar \(\lambda\) de otra manera:

\[ \hat{\lambda}(x) = max(L(\lambda | x_1, x_2, ..., x_{155}))\ ;\ \ \ \hat{\lambda}(x) = \frac{\delta}{\delta \lambda} L(\lambda | x_1, x_2, ..., x_{155}) = 0\] \[ln f(x_1, \ldots, x_{155} | \lambda) = -n · \lambda + (ln\lambda) \sum x_i - ln(\prod x_i!)\]

\[\frac{d(ln f)}{\lambda} = -n + \frac{\sum x_i}{\lambda} = 0\] \[\hat{\lambda} = \frac{\sum x_i}{n} = \bar{x}\] Comprobamos:

mean(hijos)
## [1] 1.825806

De nuevo, \(\hat{\lambda} = 1.825806\)

Estimador Bayesiano

Para el modelo Poisson existe una distribución a priori conjugada, la densidad Gamma sobre \(\lambda\), de parámetros \(a\) y \(b\).

\[f(\lambda) =\frac{b^a}{\Gamma(a)} \lambda^{a-1} exp \left( -b \lambda\right), 0 < \lambda < \infty\]

Supón que usamos una distribución inicial sobre \(\lambda \sim Ga(2,1)\). Calcula la distribución a posteriori. Dibuja con ayuda de R ambas funciones (código muy parecido al que usamos para el modelo Beta-Binomial), pero ahora hay que usar la distribución adecuada como distribución a posteriori, con los parámetros que tocan.

Proporciona un estimador puntual para \(\lambda\) usando esta distribución a posteriori.

Desarrollo del ejercicio

Histograma y función de densidad de la función a priori:

library(ggplot2)
library(magrittr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# shape = a , rate = b
gamma_values <- data.frame(rgamma(155, shape = 2, rate = 1))
colnames(gamma_values) <- c('values')

gamma_values %>%
  ggplot(aes(x=values)) +
  geom_histogram(aes(fill=..count..), binwidth = 1, col='black')

gamma_values %>%
  ggplot(aes(x=values)) +
  geom_density(alpha=.2, fill="#FF6666")

Una vez tenemos la función de distribución a priori para \(\lambda\), debemos escribir la función de verosimilitud y luego emplear el teorema de Bayes:

\[ L(x_1, x_2, ..., x_155|\lambda) = f(x_1|\lambda)·f(x_2|\lambda)·...·f(x_{155}|\lambda) \] \[ f(\lambda|x) = \frac{L(x|\lambda)·f(\lambda)}{\int{L(x|\lambda)·f(\lambda)·d\lambda}} \]

\[ P(\lambda|a,b) = \frac{b^a}{\gamma(a)} \lambda^{(a-1)}·e^{-b\lambda} \]

Donde \(\gamma(a)\) es \((a-1)!\) para enteros y \(\int_0^\infty x^{a-1}·e^{-x}dx\) en caso contrario.

\[ posterior \propto primaria · verosimilitud\]

\[ posterior \propto (e^{-n·\lambda}·\lambda^{\sum x_i}) · (\frac{b^a}{\gamma(a)} \lambda^{(a-1)}·e^{-b\lambda}) \]

\[ posterior \propto \lambda^{\sum{x_i+a-1}} e^{-\lambda · (n+b)} \]

Nota: \(e^{-n·\lambda}·\lambda^{\sum x_i}\) es el numerador de la función de verosimilitud de la función de distribución de Poisson (ver el desarrollo del ejercicio del apartado de MLE). No tenemos en cuenta el denominador porque con tener en cuenta los términos en los que aparece \(\lambda\) es suficiente.

Si desarrollamos la expresión \(\lambda^{\sum{x_i+a-1}} e^{-\lambda · (n+b)}\) veremos que se parece mucho a una distribución Gamma con parámetros \(a^* = \sum x_i + a\) y \(b^* = n + b\). Por lo tanto Gamma (que es la distribución que elegimos a priori) es una distribución conjugada.

\[ posteriori = Gamma(a^*,\ b^*) = Gamma(\sum x_i + a,\ n + b) \]

lambda_seq = seq(0, 10, length=10000)

prior_gamma_values <- data.frame(dgamma(lambda_seq, 2, 1))
colnames(prior_gamma_values) <- c('values')

post_gamma_values <- data.frame(dgamma(lambda_seq, (sum(hijos) + 2), (length(hijos) + 1)))
colnames(post_gamma_values) <- c('values')


ggplot() +
geom_line(data = prior_gamma_values, aes(x=lambda_seq, y=values)) +
geom_line(data = post_gamma_values, aes(x=lambda_seq, y=values, color='red')) +
ylim(0, 4) +
ylab('Densidad') +
xlab('Lambda') +
ggtitle('priori (negro) vs posteriori (rojo)')

gamma_post <- function(lambda) {
  return(dgamma(lambda, (sum(hijos) + 2), (length(hijos) + 1)))
}

summary(post_gamma_values)
##      values       
##  Min.   :0.00000  
##  1st Qu.:0.00000  
##  Median :0.00000  
##  Mean   :0.09999  
##  3rd Qu.:0.00000  
##  Max.   :3.69187
optimize(f = gamma_post, 
         interval = 0:10, 
         lower = 0, 
         upper = 10, 
         maximum = TRUE
         )
## $maximum
## [1] 1.820513
## 
## $objective
## [1] 3.691884

El valor de \(\lambda\) con el que más densidad se consigue (la moda de la función a posteriori) es 1.820513.

Ahora vamos a intentar calcular el error del estimador. Para ello nos vamos a basar en el método EXACT (http://onbiostatistics.blogspot.com.es/2014/03/computing-confidence-interval-for.html) que es mejor cuando las muestras son pequeñas.

number_of_events <- 1.820513 * length(hijos)

lb_exact_method <- qchisq((1-0.95)/2, 2*number_of_events)/2 
ub_exact_method <- qchisq(1-(1-0.95)/2, 2*(number_of_events+1))/2

upper_bound <- ub_exact_method / length(hijos)
lower_bound <- lb_exact_method / length(hijos)

lower_bound
## [1] 1.614276
upper_bound
## [1] 2.045796

Con un nivel de confianza del 95%, el intervalo de confianza para el \(\lambda\) estimado es \([1.614276, 2.045796]\).

Observaciones y conclusiones

Realizando el cálculo del parámetro \(\lambda\) de estas dos maneras se consigue el mismo resultado.