1. Análisis exploratorio

Los escenarios propuestos son los siguientes:

Y las distribuciones propuestas son las siguientes:

Elegimos el conjunto de datos data4_1.csv, evaluando la muestra se obtiene lo siguiente:

Estadistica Valor
Minimo 24.677791
Cuartil 1 30.975596
Mediana 33.196428
Media 33.370774
Cuartil 3 35.658315
Máximo 45.488530
Desviación Estandar 3.288815
Varianza 10.816302

Como se puede ver en la tabla, los datos se encuentran en un rango de tiempo de \(24.68\) a \(45.48\) unidades temporales. Revisando uno por uno los escenarios, intentemos ver cual hace más sentido, es decir, cuál de los escenarios puede ser descrito por la muestra.

Ahora para determinar a que distribución puede pertenecer la muestra, nos basaremos en las propiedades de las distribuciones. Se puede ver el escenario de los autos como un proceso de conteo hasta que un fenómeno ocurre n veces, ya que se cuenta el tiempo en el que pasan 100 automoviles por un punto. El conteo anterior se puede ver como un Proceso Poisson.

Recordemos que la distribución Gamma permite modelar el tiempo transcurrido hasta tener n ocurrencias de un evento generado por un proceso Poisson. Por lo tanto elegimos la distribución Gamma como la candidata a la cuál puede pertenecer la muestra.

2. Demostración estadística

Primero estimaremos los parametros de la distribución Gamma, estos son \(\beta\) (forma) y \(\lambda\) (escala), las estimaciones de los parametros son hechas siguiendo el método de máxima verosimilitud:

Parametro Valor
beta 103.336328
lambda 3.096612

A continuación usaremos 2 pruebas de bondad de ajuste, estas son la prueba Anderson-Darling y Kolmogorov-Smirnov, dichas pruebas consisten en contrastar las hipótesis:

\[H_o:\text{La muestra sigue una distribución }~Gamma(\beta,\lambda) \\ vs \\H_a:\text{La muestra no sigue una distribución }~Gamma(\beta,\lambda)\] Donde se rechaza \(H_o\) cuando el \(p-value\) de la prueba es menor al nivel de significancia \(\alpha\).

Haciendo las pruebas sobre la distribución teórica \(Gamma(\beta,\lambda)\) con \(\beta=103.34\) y \(\lambda=3.1\), para un nivel de significancia \(\alpha=0.05\) se obtiene lo siguiente:
Prueba p_value
Kolmogorov-Smirnov 0.8611279
Anderson-Darling 0.9116776

Tanto para la prueba K-S como para la prueba A-D, \(p-value>0.05\) lo que indica que no se rechaza la hipótesis nula, es decir que la muestra se ajusta a una distribución \(Gamma(\beta,\lambda)\) con \(\beta=103.34\) y \(\lambda=3.1\).

Interpretación:
\(\beta\): Es la máxima intensidad de probabilidad, este valor representa la cantidad máxima de observaciones que se esperan, en este caso \(\beta=103.34\), es decir que se espera tener a lo más \(103\) observaciones.
\(\beta\) se aproxima bastante a \(100\) que es la cantidad de coches a partir de la cual se desea estimar los tiempos de ocurrencia.

\(\lambda\): Describe el tiempo promedio entre ocurrencias del suceso, en este caso \(\lambda=3.1\) se interpreta como que en promedio, 2 automoviles pasan con \(3.1\) segundos de diferencia.

Lo anterior hace sentido ya que \(\mathbb{E}[X]=\frac{\beta}{\lambda}=\frac{103.34}{3.1}=33.3\), es decir que en promedio se espera que los \(100\) coches pasen en \(33.3\) segundos, lo que es igual a decir que pase un coche cada 3 segundos.

3. Representación visual

A continuación se muestran algunas gráficas que permitirán ver la relación de la muestra y la distribución teórica.

Todo lo anterior refuerza y complementa la evidencia estadistica, por lo que se concluye sin dudas que la muestra sigue una distribución \(Gamma\).

4. Comparación

Comparemos los valores estadisticos muestrales obtenidos al inicio, con los valores estadisticos poblacionales teóricos.

Estadistico Valor muestral Valor teórico
Minimo 24.677791 0.0019429
Cuartil 1 30.975596 31.1018543
Mediana 33.196428 33.2631916
Media 33.370774 33.3707742
Cuartil 3 35.658315 35.5224422
Máximo 45.488530 0.0004007
Desviación Estandar 3.288815 3.2827648
Varianza 10.816302 10.7765448

Como se puede observar los valores teóricos se asemejan bastante a los valores muestrales, más en el caso de la media y la desviación estandar donde la precisión es de 2 decimales, en el caso de los demás valores la diferencia entre el teórico y el muestral no supera a \(0.2\).

Para el mínimo y el máximo, los valores que se colocaron en la columna Valor teórico son: \[\mathbb{P}[x\leq minimo~muestral]~~ \text{ y }~~ \mathbb{P}[x\geq maximo~muestral]\] respectivamente, y como se puede observar estas probabilidades son muy pequeñas, lo que indica que las colas de la distribución se ajustan bastante bien con los extremos muestrales, que todo valor que esté fuera de este intervalo (min,max) tendrá una probabilidad insignificante y poco representativa.

5. Función de supervivencia, de riesgo y de riesgo acumulado

Calcula la función de riesgo y función de riesgo acumulado como si se tratara del caso discreto y compararlos con los teóricos.

Se observa que la función de supervivencia muestral parece ajustarse bien a la teórica aunque en la parte central de forma casi imperceptible no se acerca tanto, en cambio las funciones \(h(t)\) y \(H(t)\) parecen distanciarse cada vez más de la función teórica, esto debido a que son funciones sensibles a los cambios de \(S(t)\), las pequeñas diferencias entre la \(S(t)\) muestral y la teórica generan diferencias más significativas en las funciones de riesgo y por lo tanto de riesgo acumulado.

6. Extra

Usando las observaciones con censura por la izquierda y por la derecha, procedemos a calcular los estimadores máximo verosimil para ambos casos censurados:

Parametro Muestra Normal Censura Derecha Censura Izquierda
beta 103.336328 94.977041 94.197296
lambda 3.096612 2.797496 2.864186

No hay mucho que añadir, los parametros estimados no difieren mucho entre sí ni de los parametros de la muestra con las observaciones completas.

A continuación se grafican las distribuciones teóricas de las muestras a partir de los parametros estimados.

Debido a que los parametros no difieren entre si, en consecuencia las gráficas de las densidades teóricas tampoco difieren mucho entre sí. Para los datos censurados por la izquierda la gráfica se desplaza un poco hacia la izquierda y de la misma manera para los datos censurados por la derecha, la gráfica se desplaza hacia la derecha.

Anexo

Código usado

# 1. ----

#Importamos los datos
datos4<-read.csv("data4_1.csv", header = FALSE)
datos4<-data.frame(datos4)
d<-datos4$V1 #<- vector que contiene los datos de la muestra

library(knitr)
library(kableExtra)

#generamos una tabla con los valores estadisticos más significativos de la muestra
est<-c("Minimo","Cuartil 1","Mediana","Media","Cuartil 3","Máximo","Desviación Estandar","Varianza")
val<-c(min(d),quantile(d, 0.25),median(d),mean(d),quantile(d, 0.75),max(d),sd(d),var(d))

#construimos una tabla para graficar los valores
kable(data.frame(Estadistica=est,Valor=val)) %>% 
    kable_styling(full_width = F)

#  2.----

library(univariateML)
#Estimamos los parametros de la distribucion por maxima verosimilitud
gam<-mlgamma(datos4$V1)

#almacenamos los parametros en un vector
param<-c(gam[[1]],gam[[2]])
par<-c("beta","lambda")
#desplegamos la tabla de valores
kable(data.frame(Parametro=par,Valor=param)) %>% 
    kable_styling(full_width = F)

library(goftest)
#realizamos las pruebas estadisticas Kolmogorov-Smirnov y Anderson-Darling
Ks<- ks.test(d, "pgamma", shape=gam[[1]],rate=gam[[2]])
Ad<- ad.test(d, "pgamma", shape=gam[[1]],rate=gam[[2]])

prueba<-c("Kolmogorov-Smirnov","Anderson-Darling")
pval<-c(Ks$p.value,Ad$p.value)
kable(data.frame(Prueba=prueba,p_value=pval)) %>% 
    kable_styling(full_width = F)

# 3.----
library(ggplot2)
require(gridExtra)
library(qqplotr)
#histograma + distribucion
p1<-ggplot(datos4, aes(x=V1))+   
  geom_histogram(binwidth=1, aes(y =..density..), 
                 col='black',fill='steelblue', alpha=0.9)+
  stat_function(fun = function(.x){dml(x = .x, obj = mlgamma(datos4$V1))},
                aes(color = "gamma"), size = 1) +
  scale_color_manual(breaks = "gamma",
                     values = c("gamma" = "red")) +
  labs(color = "Distribución") +
  theme_bw() +
  theme(legend.position = "bottom")+
  labs(title="Muestra vs Distribución Teórica", x="Valor", y="Frecuencia")+
  theme(plot.title = element_text(hjust = 0.5))

# QQ-Plot
p2<-ggplot(datos4, aes(sample = V1)) +
  stat_qq_point(distribution = "gamma", dparams = c(gam[[1]],gam[[2]]), col="steelblue") +
  stat_qq_line(distribution = "gamma", dparams = c(gam[[1]],gam[[2]]))+
  labs(title="QQ-Plot", x="Cuantiles teóricos", y="Cuantiles muestrales")+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
#PP-Plot
p3<-ggplot(datos4, aes(sample = V1)) +
  stat_pp_point(distribution = "gamma", dparams = list(gam[[1]],gam[[2]]), col="steelblue") +
  stat_pp_line(distribution = "gamma", dparams = c(gam[[1]],gam[[2]]))+
  labs(title="PP-Plot",x="Cuantiles teóricos", y="Cuantiles muestrales")+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
#Funcion de distribucion acumulada
p4<-ggplot(datos4, aes(x = V1))+
  stat_ecdf(geom = "step", pad = FALSE)+
  stat_function(fun=pgamma,color="steelblue",args=list(gam[[1]],gam[[2]]), 
                aes(color = "gamma"), size = 1)+
  scale_color_manual(breaks = "gamma",
                     values = c("gamma" = "steelblue")) +
  labs(color = "Distribución")+
  labs(title="F(x) Teórica vs Empirica",x="Valor", y="Probabilidad")+
  theme_bw() +
  theme(legend.position = "bottom")+
  theme(plot.title = element_text(hjust = 0.5))

#desplegamos las 4 graficas juntas
grid.arrange(p1,p2,p3,p4, ncol=2, nrow=2)

# 4.----
#calculamos los valores estadisticos reales de la distribucion
min<-pgamma(24.677791, gam[[1]], gam[[2]], lower.tail = T, log.p = F)   
q1<-qgamma(0.25, gam[[1]],gam[[2]], lower.tail = T, log.p = F)
m<-gam[[1]]/gam[[2]]
me<-qgamma(0.5, gam[[1]],gam[[2]], lower.tail = T, log.p = F)
q3<-qgamma(0.75, gam[[1]],gam[[2]], lower.tail = T, log.p = F)
max<-pgamma(45.48853, gam[[1]], gam[[2]], lower.tail = F, log.p = F)
sd<-sqrt(gam[[1]]/gam[[2]]^2)
var<-gam[[1]]/gam[[2]]^2

est<-c("Minimo","Cuartil 1","Mediana","Media","Cuartil 3","Máximo","Desviación Estandar","Varianza")
val<-c(min(d),quantile(d, 0.25),median(d),mean(d),quantile(d, 0.75),max(d),sd(d),var(d))
rval<-c(min,q1,me,m,q3,max,sd,var)

comparacion<-data.frame(est, val, rval)
colnames(comparacion)<- c("Estadistico","Valor muestral","Valor teórico")
# tabla comparativa
kable(comparacion) %>% 
    kable_styling(full_width = F)

# 5. ----
library(flexsurv)
#leemos el archivo que contiene la funcion F(x) de la muestra
datos<-read.csv("data4_2.csv", header = FALSE)

l<-length(datos$V1)
t<-seq(1,l)

#funciones teóricas
s_t<-1-pgamma(t, gam[[1]],gam[[2]])
h_t<-hgamma(t,gam[[1]],gam[[2]])
H_t<-Hgamma(t,gam[[1]],gam[[2]])

#obtenemos las funciones muestrales
s_m<-1-(datos[[1]])
library(dplyr)
h_m= 1-(lead(s_m,offset = 1, defaultValue = 0)/s_m)
H_m=cumsum(h_m)

#graficamos las funciones S,h y H
p1<-ggplot() + 
  geom_line(aes(x=t,y=s_m, color="muestral"), size=1) + 
  geom_line(aes(x=t,y=s_t, color="teorica"),size=1)+
  scale_color_manual(breaks = c("muestral", "teorica"),
                     values = c("muestral" = "steelblue", "teorica" = "red")) +
  theme_bw() +
  labs(title="Funcion de Supervivencia", x="Tiempo", y="S(t)",color = "Distribución")+
  theme(plot.title = element_text(hjust = 0.5))  
  
  
p2<-ggplot() + 
  geom_line(aes(x=t,y=h_m, color="muestral"), size=1) + 
  geom_line(aes(x=t,y=h_t, color="teorica"),size=1)+
  scale_color_manual(breaks = c("muestral", "teorica"),
                     values = c("muestral" = "steelblue", "teorica" = "red")) +
  theme_bw() +
  labs(title="Funcion de Riesgo", x="Tiempo", y="h(t)",color = "Distribución")+
  theme(plot.title = element_text(hjust = 0.5))

p3<-ggplot() + 
  geom_line(aes(x=t,y=H_m, color="muestral"), size=1) + 
  geom_line(aes(x=t,y=H_t, color="teorica"),size=1)+
  scale_color_manual(breaks = c("muestral", "teorica"),
                     values = c("muestral" = "steelblue", "teorica" = "red")) +
  theme_bw() +
  labs(title="Funcion de Riesgo Acumulado", x="Tiempo", y="H(t)",color = "Distribución")+
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(p1,p2,p3, ncol=2, nrow=2)

# 6.----

library(survival)
library(dplyr)
library(fitdistrplus)
#cargamos los archivos que contienen las censuras
censura1<-read.csv("data4_left_censored.csv", header = FALSE)
censura2<-read.csv("data4_right_censored.csv", header = FALSE)
censura1<-data.frame(censura1)
censura2<-data.frame(censura2)

censura_derecha<-censura2$V1
censura_izquierda<-censura1$V1

#d es la muestra
#generamos un objeto de tipo Surv, con censura por la derecha y por la izquierda
datos_censurados_derecha <- Surv(d,censura_derecha)
datos_censurados_izquierda <- Surv(d, censura_izquierda, type = "left")

dc_derecha<- tibble(left = datos_censurados_derecha[,"time"]) %>% 
  mutate(right = if_else(datos_censurados_derecha[,"status"] == 0, left, NA_real_))
#obtenemos los parametros ajustados a la censura por maxima verosimilitud
ajuste <- fitdistcens(as.data.frame(dc_derecha), "gamma")

dc_izquierda <- tibble(right = datos_censurados_izquierda[,"time"]) %>% 
  mutate(left = if_else(datos_censurados_izquierda[,"status"] == 0, right, NA_real_))
#obtenemos los parametros ajustados a la censura por maxima verosimilitud
ajuste2 <- fitdistcens(as.data.frame(dc_izquierda), "gamma")

param<-c("beta","lambda")
param1<-c(gam[[1]],gam[[2]])
param2<-c(ajuste$estimate[[1]],ajuste$estimate[[2]])
param3<-c(ajuste2$estimate[[1]],ajuste2$estimate[[2]])
parametros<-data.frame(param,param1,param2,param3)
colnames(parametros)<- c("Parametro","Muestra Normal","Censura Derecha", "Censura Izquierda")
kable(parametros) %>% 
  kable_styling(full_width = F)

# grafica para comparar las 3 distribuciones
ggplot(datos4, aes(x=V1))+   
  stat_function(fun = ~dgamma(.x, shape= gam[[1]], 
                              rate = gam[[2]]), aes(color = "Sin censura"), size = 1)+
  stat_function(fun = ~dgamma(.x, shape= ajuste$estimate[[1]], rate = ajuste$estimate[[2]]), 
                aes(color = "Censura derecha"), size = 1)+
  stat_function(fun = ~dgamma(.x, shape= ajuste2$estimate[[1]], rate = ajuste2$estimate[[2]]), 
                aes(color = "Censura izquierda"), size = 1)+
  scale_color_manual(breaks = c("Sin censura", "Censura derecha","Censura izquierda"),
                     values = c("Sin censura" = "steelblue", "Censura derecha" = "red",
                                "Censura izquierda"="green")) +
  labs(color = "Distribución") +
  theme_bw() +
  theme(legend.position = "bottom")+
  labs(title="Distribuciones con Censura", x="Valor", y="Probabilidad")+
  theme(plot.title = element_text(hjust = 0.5))