Julian Andres Ospina Cuesta

Maestría en Ciencia de Datos - Universidad Javeriana Cali
Cali, Colombia
Código: 8984933 +57 3043912471
library(ggplot2)
library(dplyr)
library(knitr)
library(gridExtra)

Problema 1: Estimación del valor de π

En el siguiente problema se desea estimar el valor \(π\) con una simulación que consta de identificar los puntos dentro del circulo que esta a su vez dentro de un cuadrado de área igual a 1. se deben elegir \(n\) puntos distribuidos en el cuadrado y se sabe que \(f( punto dentro del circulo) = π/4\).

Para desarrollar este problema vamos a seguir los pasos que sugieré el planteamiento con \(n\) valores generados:

#Guardar mi respuesta en cada corrida
df <- data.frame(n = character(0), valor_pi = numeric(0), proporcion = numeric(0))

#configuración de notación cientifica
options(scipen = 10)

#Nota: con fines de replicar el ejercicio, se fija la semilla generadora de números
set.seed(100)

1. Calculos con diversos tamaños de datos generados

n=1000

1.1. Generación de datos \(X∼unif(x,a=0,b=1)\) tanto ‘X’ como para ‘Y’

Cada valor en \(X_i\) y \(Y_i\) es una coordena de un punto dentro de mi cuadrado

n = 1000

X_muestral=runif(n, min=0, max=1) 
Y_muestral=runif(n, min=0, max=1)

1.2. Calculo de distancia al centro del cuadrado \((X_{centro}=0.5, y_{centro}=0.5)\)

Se pide emplear la función de distancia \(f(x,y) = (X_i − 0.5)^2 + (Y_i − 0.5)^2\), e identificar que sea menor que 0.25

df_muestral <- data.frame(x = X_muestral, y = Y_muestral)
df_muestral$distancia <- apply(df_muestral, 1, function(row) (row["x"] - 0.5)^2 + (row["y"] - 0.5)^2)
df_muestral$en_circulo <- ifelse(df_muestral$distancia <= 0.25, 1, 0)

head(df_muestral)
##            x         y  distancia en_circulo
## 1 0.30776611 0.0740483 0.21838872          1
## 2 0.25767250 0.1118664 0.20937034          1
## 3 0.55232243 0.6239440 0.01809975          1
## 4 0.05638315 0.6710818 0.22606489          1
## 5 0.46854928 0.3658942 0.01897351          1
## 6 0.48377074 0.1831814 0.10063742          1

1.3.Identifcar el número de puntos adentro y estimar el valor de π

Veamos los siguientes enunciados:

  • El \(Area_{circulo} = π/4\), eso significa que sabiendo el área de circulo podemos calcular \(π\)
  • El circulo esta dentro de un cuadrado de área 1, es decir que si sacamos proporciones del cuadrado, podemos estimar de manera proporcional el área del circulo.
  • La proporción de puntos dentro del circulo sobre el total de puntos, me dará la proporcion del área que requeriamos para estimar el área del circulo, y posteriormente me ayudará con la estimación \(π\)
proporcion_circulo <- sum(df_muestral$en_circulo)/n
area_circulo = 1 * proporcion_circulo #Recordemos que el área del cuadrado es 1
valor_pi_estimado =  area_circulo*4

df <- rbind(df, data.frame(n=as.character(n), valor_pi =valor_pi_estimado, proporcion =proporcion_circulo)) 

cat('Proporción_circulo:',proporcion_circulo)
## Proporción_circulo: 0.775
cat('Valor_pi:',valor_pi_estimado)
## Valor_pi: 3.1

n=10000 y n= 100000

Definir funcion para n cantidad

estimacion_valor_pi<-function(n){
  #función empleada para el calculo de estimados del valor pi, a partir de la generación
  #de **n** tamaño definido.
  #**NOTA**: Esta función es la autoamtización de calculos de la pestaña n=1000
  #
  #Parametros:
  #@n(integer): cantidad de números a generar
  #
  #Salida
  #@return(dataframe): data frame con un dataframe resumen
  
  # 1. Generación de datos
  X_muestral=runif(n, min=0, max=1) 
  Y_muestral=runif(n, min=0, max=1) 
  
  # 2. Calculo de distancia al centro del cuadrado
  df_muestral <- data.frame(x = X_muestral, y = Y_muestral)
  df_muestral$distancia <- apply(df_muestral, 1, function(row) (row["x"] - 0.5)^2 + (row["y"] - 0.5)^2)
  df_muestral$en_circulo <- ifelse(df_muestral$distancia <= 0.25, 1, 0)
  
  # 3. Identifcar el número de puntos adentro y estimar el valor de π
  proporcion_circulo <- sum(df_muestral$en_circulo)/n
  area_circulo = 1 * proporcion_circulo #Recordemos que el área del cuadrado es 1
  valor_pi_estimado =  area_circulo*4

  return (data.frame(n=as.character(n), valor_pi =valor_pi_estimado, proporcion =proporcion_circulo)) 
}

Calculando estimador para los n definidos

n_requeridos = c(10000,100000)
for (n_muestral in n_requeridos) {
  df <- rbind(df,estimacion_valor_pi(n_muestral))
}

2. Conclusión del estimador del valor π

Veamos los datos recolectados

##        n valor_pi proporcion
## 1   1000  3.10000    0.77500
## 2  10000  3.14400    0.78600
## 3 100000  3.14204    0.78551

Se observa una clara relación entre el tamaño de n y la estimación de \(π\), es decir, Entre mas se aumenta el tamaño de la muestra , mayor es el acercamiento a la proporción del estimada del área de circulo a partir del cuadrado y por consiguiente mayor es la aproximación al valor \(π\).

Problema 2: Propiedades de los estimadores

Se nos piden realizar 4 estimadores para una posterior validación de propiedades, estos estimadores se obtienen con la generación de 4 valores que siguen una distribucón \(X \sim \text{Exp}(\lambda = 1/θ_{supuesto})\)

En el ejercicio se nos pide soponer un parametro \(θ\) y una muestra de maximo 1000 estimadores

#Guardar mi respuesta en cada corrida
Theta_supuesto = 10
la_matrix <-matrix(data=rexp(1000*4, rate=1/Theta_supuesto), nrow=1000, byrow=TRUE)
df_datos <- as.data.frame(la_matrix)
colnames(df_datos) <- c("X1", "X2", "X3","X4")

#Para los futuros calculos, tengan una descripción
conjunto_nombres = c('theta_1','theta_2','theta_3','theta_4')

A continuación observamos los estimadores que tenemos durante este ejercicio:

\[\hat{θ_1}=\frac{X_1+X_2}{6}+\frac{X_3+X_4}{3}\]

\[\hat{θ_2}=\frac{X_1+2∗X_2+3∗X_3+4∗X_4}{5}\]

\[\hat{θ_3}=\frac{X_1+X_2+X_3+X_4}{4}\]

\[\hat{θ_4}=\frac{min(X_1,X_2,X_3,X_4)+max(X_1,X_2,X_3,X_4)}{2}\]

1. Calculo de los diferentes \(\hat{θ_{1,2,3,4}}\)

theta_1 <- function(x){
  return( (x[1] + x[2]) / 6 + (x[3] + x[4]) / 3 )
}

theta_2 <- function(x){
  return( (x[1] + 2 * x[2] + 3 * x[3] + 4 * x[4] ) / 5 )
}

theta_3 <- function(x){
  return( (x[1] + x[2] + x[3] + x[4]) / 4 )
}

theta_4 <- function(x){
  return( (min(x) + max(x)) / 2 )
}

#Creamos mis estimadores en mi muestra de datos
df_datos <- df_datos %>%
  mutate(theta_1 = apply(df_datos, 1, theta_1))
df_datos <- df_datos %>%
  mutate(theta_2 = apply(df_datos, 1, theta_2))
df_datos <- df_datos %>%
  mutate(theta_3 = apply(df_datos, 1, theta_3))
df_datos <- df_datos %>%
  mutate(theta_4 = apply(df_datos, 1, theta_4))

2. Propiedades de los estimadores

A continuación planteamos las funciones que nos permitiran calcular cada uno de las proppiedades

Insesgado

Para calcular el sesgado de nuestros datos, podemos emplear la propiedad

\[E(\hat{ \Theta}) = \Theta\]

ó

\[ \Theta - E(\hat{ \Theta})= 0\]

insesgado <- function(df_seccion){
  #Funcion empleada para calcular en las estimaciones el insesgado (que viene por columna en el df (1-4)),
  #y devolver el valor esperado por los 4 estimadores en un conjunto. 
  #Nota: La idea es que nos de muy cerca de 0
  #
  #Parametros:
  #@df_seccion(dataframe): el dataframe con los datos previamente seccionados
  #(Es decir, separados por el k requerido para el analisis, pero con los 4 estimadores por columna)
  #
  #Salida
  #@return(conjunto): data frame con un dataframe resumen

  estimaciones<-c()
  for (theta_calculado in conjunto_nombres) {
    estimaciones <-append(estimaciones,( Theta_supuesto - mean(df_seccion[[theta_calculado]])))
  }
  return(estimaciones)
}

Eficiencia

Para calcular la eficiencia debemos apoyarnos en el concepto de eficiencia relativa \[eff(\hat{ \theta_1},\hat{ \theta_2}) = \frac{Var(\hat{ \theta_2})}{Var(\hat{ \theta_1})}\] Es intereseante ver que siempre comparo dos variaciones, pero, ¿y si comparamos la variacion esperada de mi población?, la cual la conocemos gracias a mi distribución de datos conocida, la exponencial en donde sabemos:

\[\frac{1}{\lambda} = θ_{supuesto} = \mu\] y sabemos que tambien que se cumple en una distribución exponencial

\[ \mu^2 = \sigma^2 \] y podria ser interesante mirar mi eficiencia frente variación de mi población

eficiencia <- function(df_seccion){
  #Funcion empleada para calcular en las estimaciones la eficiencia (que viene por columna en el df (1-4)),
  #y devolver el valor de la eficiencia relativa por los 4 estimadores en un conjunto
  #Nota: La idea es que nos de un número mayor de 1
  #
  #Parametros:
  #@df_seccion(dataframe): el dataframe con los datos previamente seccionados
  #(Es decir, separados por el k requerido para el analisis, pero con los 4 estimadores por columna)
  #
  #Salida
  #@return(conjunto): data frame con un dataframe resumen

  estimaciones<-c()
  for (theta_calculado in conjunto_nombres) {
    estimaciones <-append(estimaciones,Theta_supuesto^2 /var(df_seccion[[theta_calculado]]))
  }
  return(estimaciones)
}

Consistencia

Para calcular la concistencia nos ayudaremos del Error Cuadratico Medio

\[EMC(\hat{ \theta})= Var(\hat{ \theta}) + sesgo(\hat{ \theta},\theta)^2\]

consistencia <- function(df_seccion){
  #Funcion empleada para calcular en las estimaciones la consistencia (que viene por columna en el df(1-4)),
  #y devolver el valor EMC(Error cuadrático Medio) por los 4 estimadores en un conjunto
  #
  #Parametros:
  #@df_seccion(dataframe): el dataframe con los datos previamente seccionados
  #(Es decir, separados por el k requerido para el analisis, pero con los 4 estimadores por columna)
  #
  #Salida
  #@return(conjunto): data frame con un dataframe resumen

  estimaciones<-c()
  for (theta_calculado in conjunto_nombres) {
    estimaciones <-append(estimaciones,var(df_seccion[[theta_calculado]])+( mean(df_seccion[[theta_calculado]])- Theta_supuesto)^2)
  }
  return(estimaciones)
}

3. Generar las n separaciones solicitadas

Empleamos nuestro dataframe de datos y separamos en las muestras solicitadas

n_separaciones = c(20, 50, 100,1000)

dfs_particion <- list()
for (n_separacion in n_separaciones) {
  df_temporal<-df_datos[1:n_separacion,conjunto_nombres]
  dfs_particion <-append(dfs_particion,list(df_temporal))
}

4. Calculemos las propiedades de los estimadores

dfs_analisis = list()
for (df_part in dfs_particion) {
  conj_insesgado <- insesgado(df_part)
  conj_eficiencia <- eficiencia(df_part)
  conj_consistencia <- consistencia(df_part)
  df <- data.frame(thetas = conjunto_nombres,
                   insesgado = conj_insesgado,
                   eficiencia = conj_eficiencia,
                   consistencia =conj_consistencia )
  dfs_analisis <-append(dfs_analisis,list(df))
}

5. Revision de resultados

n = 20

df = dfs_particion[[1]]
boxplot(df, main = "Diagrama de caja de bigotes de Theta's'", xlab = "Theta", ylab = "Valores")

kable(do.call(data.frame, dfs_analisis[[1]]),
      format = "markdown",  
      caption = "**Tabla n°1 Resumen de estimados con n=20**",
      align = "c", escape = FALSE,
      row.names = FALSE,
      booktabs = TRUE)
Tabla n°1 Resumen de estimados con n=20
thetas insesgado eficiencia consistencia
theta_1 -0.3024511 4.427134 22.67946
theta_2 -11.2598883 0.941090 233.04485
theta_3 0.1846227 6.029423 16.61942
theta_4 -2.9160952 3.152022 40.22928

Se observa que el \(\hat{θ_3}\) es la estimación que mas insesgada (más se acerca al cero) y menor cifra de EMC,indicador que nos indica la consistencia. tambien debemos añadir que su eficiencia es la mas alta, eso se debe a que su variación es menor que el los demas. Tambien se menciona que el \(\hat{θ_2}\) no es un buen estimador, vemos mucha variación que se regleja en los peores resultados dentro de mis estimadores.

n = 50

df = dfs_particion[[2]]
boxplot(df, main = "Diagrama de caja de bigotes de Theta's'", xlab = "Theta", ylab = "Valores")

kable(do.call(data.frame, dfs_analisis[[2]]),
      format = "markdown",  
      caption = "**Tabla n°1 Resumen de estimados con n=50**",
      align = "c", escape = FALSE,
      row.names = FALSE,
      booktabs = TRUE)
Tabla n°1 Resumen de estimados con n=50
thetas insesgado eficiencia consistencia
theta_1 -0.5347040 4.0675727 24.87060
theta_2 -11.7763907 0.8211207 260.46815
theta_3 -0.2608704 4.7109005 21.29542
theta_4 -3.4711825 2.5817111 50.78311

Se observa que el \(\hat{θ_3}\) vuelve a destacar, pero vemos que ya empieza el boxplot a mostrar datos atipicos y que en general las colas de todos los estimadores son muy largas.

n = 100

df = dfs_particion[[3]]
boxplot(df, main = "Diagrama de caja de bigotes de Theta's'", xlab = "Theta", ylab = "Valores")

kable(do.call(data.frame, dfs_analisis[[3]]),
      format = "markdown",  
      caption = "**Tabla n°1 Resumen de estimados con n=100**",
      align = "c", escape = FALSE,
      row.names = FALSE,
      booktabs = TRUE)
Tabla n°1 Resumen de estimados con n=100
thetas insesgado eficiencia consistencia
theta_1 0.0747816 5.046831 19.82001
theta_2 -10.4202483 1.042935 204.46481
theta_3 0.1047914 5.392149 18.55646
theta_4 -2.8329314 3.062494 40.67863

Se observa que el \(\hat{θ_3}\) sigue con al tendencia de ser el mejor estimador, pero ahora vemos que tenemos tendencia de ir bajando variacion en mis thetas si comparamos con el tamaño de n =50

n = 1000

df = dfs_particion[[4]]
boxplot(df, main = "Diagrama de caja de bigotes de Theta's'", xlab = "Theta", ylab = "Valores")

kable(do.call(data.frame, dfs_analisis[[4]]),
      format = "markdown",  
      caption = "**Tabla n°1 Resumen de estimados con n=1000**",
      align = "c", escape = FALSE,
      row.names = FALSE,
      booktabs = TRUE)
Tabla n°1 Resumen de estimados con n=1000
thetas insesgado eficiencia consistencia
theta_1 -0.3582129 3.3923891 29.60607
theta_2 -10.7043033 0.7828883 242.31425
theta_3 -0.3341044 3.7172455 27.01327
theta_4 -3.2054664 2.0910343 58.09824

Se observa que el \(\hat{θ_3}\) es en definitiva nuestro estimador estrella, siempre se aproximo a nuestro \(θ_{supuesto}=10\) y que vemos como las colas se han acortado en todos los thetas, apareciendo mayor cantidad de datos atipicos, y hablandonos sobre lo consistentes que se van volviendo nuestros estimadores al incrementar el \(n\) en nuestros muestreos.

Problema 3: Teorema Central

Se nos pide analizar una población de 1000 plantas con proporción de enfermas cambiantes, para estudiar la estimación proporción. mirando mediana, varianza y normalidad al seleccionar diversos tamaños de muestras.

Vamos a iniciar creando nuestra función de creación de población con base en una proporción de enfermos

1. Creando la población

N_poblacion = 1000

obtener_poblacion <- function(proporcion_enferma){
  #Esta función devuelve una proporción de 1 y 0 que simboliza pacientes enfermos (1) y sanos (0) a partir de una proporción dada.
  #Nota: No se emplea rbinom(N_poblacion,1,proporcion_enferma) como sugiere el problema en las notas, por sugerencia dada en clase del 6 sept
  #
  #Parametros:
  #@proporcion_enferma(float): Es el % de enfermos que tiene la población que quiero obtener
  #
  #Salida
  #@return(conjunto): regresa un conjunto de datos que representan mi población
  
  cantidad_enfermos = round(N_poblacion*proporcion_enferma)
  return(append(replicate(cantidad_enfermos, 1), replicate(N_poblacion-cantidad_enfermos, 0)))
}

2. Generando la muestra con la proporción \(\hat{p}\)

obtener_muestra<- function(proporcion_enferma,n_muestra){
  
  Poblacion <- obtener_poblacion(proporcion_enferma)
   #Nota:El enunciado no especifica si el muestreo es con reemplazo, se asume que no. dado que son espécimenes de plantas
  la_muestra <- sample(Poblacion,n_muestra,replace=TRUE)
  return (mean(la_muestra))#estimado p gorro se puede ver como media de la muestra obtenida.
}

3. Calculando indicadores al estimador en las muestras solicitadas

El ejercicio solicita hacer un analisis de la muestra frente a su proporción, luego nos pide que a dicho estimador calculemos la varianza,mediana y posteriormente de una prueba de shapiro para ver normalidad. Por lo que es conveniente crear una función que calcule los estimados dado las repeticiones del muestreo que nos soliciten.

realizar_analisis<-function(proporcion_enferma,n_muestra,repeticion_muestreo){
  
  La_muestra <- c()
  for (i in 1:repeticion_muestreo){
    La_muestra <-append(La_muestra,obtener_muestra(proporcion_enferma,n_muestra))
  }
  mediana=mean(La_muestra)
  varianza <- var(La_muestra)
  sesgo <- mediana - proporcion_enferma
  valor_p_shapiro <- shapiro.test(La_muestra)$p.value
  h0_normalidad <- valor_p_shapiro>0.05
  
  return(list(La_muestra,
    data.frame(n_muesra = n_muestra, mediana=mediana, varianza=varianza, sesgo = sesgo,
                    valor_p_shapiro= valor_p_shapiro, normalidad = h0_normalidad)))   
}

4. Analisis en diversos tamaños de muestra

Veamos como se comportan nuestra investigacion con los siguientes tamaños n = 5 10 15 20 30 50 60 100 200 500 con una repetición del ejercicio de 500 veces por cada n.

Proporción del 50%

n_solicitados = c(5,10,15,20,30,50,60,100,200,500)
repeticion_muestreo = 500
proporcion_enferma = 50/100 #50%

df <- data.frame()

#Configurar el diseño de una ventana gráfica dividida en múltiples paneles o celdas
par(mfrow=c(1,2))
for (n_muestra in n_solicitados){
  resultados = realizar_analisis(proporcion_enferma,n_muestra,repeticion_muestreo)
  hist(resultados[[1]],freq=FALSE,col="blue",main=paste("Histograma de tamaño: ",n_muestra),xlab = "Proporción")
  qqnorm(resultados[[1]],main=paste("QQ-plot de tamaño: ",n_muestra))
  qqline(resultados[[1]], col="Blue")

  df <-rbind(df,resultados[[2]])

}

#Configurar el diseño de una ventana gráfica dividida en múltiples paneles o celdas
par(mfrow=c(1,1))

kable(do.call(data.frame, df),
      format = "markdown",  
      caption = "**Tabla n°2 Resumen de estimados con 50% de enfermos**",
      align = "c", escape = FALSE,
      row.names = FALSE,
      booktabs = TRUE)
Tabla n°2 Resumen de estimados con 50% de enfermos
n_muesra mediana varianza sesgo valor_p_shapiro normalidad
5 0.5060000 0.0510261 0.0060000 0.0000000 FALSE
10 0.5050000 0.0259269 0.0050000 0.0000000 FALSE
15 0.4961333 0.0168276 -0.0038667 0.0000004 FALSE
20 0.5010000 0.0134359 0.0010000 0.0000013 FALSE
30 0.5038000 0.0090815 0.0038000 0.0001362 FALSE
50 0.5004400 0.0050235 0.0004400 0.0012233 FALSE
60 0.4965000 0.0044739 -0.0035000 0.0159990 FALSE
100 0.5010000 0.0024230 0.0010000 0.0090343 FALSE
200 0.5004800 0.0013729 0.0004800 0.1632933 TRUE
500 0.5000480 0.0004739 0.0000480 0.2917598 TRUE

Con la proporcion de 50% vemos que los histogramas tratan de ser muy simetricos, y a medida que subimos de n en la muestra, la hipotesis nula se anula y podemos concluir que los valores se comportan normales. su indicadores de tendencia central y de disperción mejoran a medidas que incrementamos n.

Proporción del 10%

n_solicitados = c(5,10,15,20,30,50,60,100,200,500)
repeticion_muestreo = 500
proporcion_enferma = 10/100 #10%

df <- data.frame()
par(mfrow=c(1,2))
for (n_muestra in n_solicitados){
  resultados = realizar_analisis(proporcion_enferma,n_muestra,repeticion_muestreo)
  hist(resultados[[1]],freq=FALSE,col="blue",main=paste("Histograma de tamaño: ",n_muestra),xlab = "Proporción")
  qqnorm(resultados[[1]],main=paste("QQ-plot de tamaño: ",n_muestra))
  qqline(resultados[[1]], col="Blue")

  df <-rbind(df,resultados[[2]])

}

par(mfrow=c(1,1))

kable(do.call(data.frame, df),
      format = "markdown",  
      caption = "**Tabla n°3 Resumen de estimados con 10% de enfermos**",
      align = "c", escape = FALSE,
      row.names = FALSE,
      booktabs = TRUE)
Tabla n°3 Resumen de estimados con 10% de enfermos
n_muesra mediana varianza sesgo valor_p_shapiro normalidad
5 0.1024000 0.0169081 0.0024000 0.0000000 FALSE
10 0.0968000 0.0084467 -0.0032000 0.0000000 FALSE
15 0.1012000 0.0054049 0.0012000 0.0000000 FALSE
20 0.1042000 0.0047619 0.0042000 0.0000000 FALSE
30 0.1043333 0.0029716 0.0043333 0.0000000 FALSE
50 0.1004400 0.0018796 0.0004400 0.0000000 FALSE
60 0.0975000 0.0013693 -0.0025000 0.0000007 FALSE
100 0.0999400 0.0008647 -0.0000600 0.0000937 FALSE
200 0.1007300 0.0004779 0.0007300 0.0003560 FALSE
500 0.1000160 0.0001723 0.0000160 0.0405607 FALSE

Con la proporcion de 10% vemos que los histogramas al inicio no son muy simetricos, pero con el tiempo casi se llega a una distribución normal segun la prueba shapiro con el mayor n en el ejercicio. su indicadores de tendencia central y de disperción mejoran a medidas que incrementamos n.

Proporción del 90%

n_solicitados = c(5,10,15,20,30,50,60,100,200,500)
repeticion_muestreo = 500
proporcion_enferma = 90/100 #90%
graficos<- list()

df <- data.frame()
par(mfrow=c(1,2))
for (n_muestra in n_solicitados){
  resultados = realizar_analisis(proporcion_enferma,n_muestra,repeticion_muestreo)
  hist(resultados[[1]],freq=FALSE,col="blue",main=paste("Histograma de tamaño: ",n_muestra),xlab = "Proporción")
  qqnorm(resultados[[1]],main=paste("QQ-plot de tamaño: ",n_muestra))
  qqline(resultados[[1]], col="Blue")

  df <-rbind(df,resultados[[2]])

}

par(mfrow=c(1,1))

kable(do.call(data.frame, df),
      format = "markdown",  
      caption = "**Tabla n°4 Resumen de estimados con 90% de enfermos**",
      align = "c", escape = FALSE,
      row.names = FALSE,
      booktabs = TRUE)
Tabla n°4 Resumen de estimados con 90% de enfermos
n_muesra mediana varianza sesgo valor_p_shapiro normalidad
5 0.9060000 0.0189619 0.0060000 0.0000000 FALSE
10 0.8948000 0.0085901 -0.0052000 0.0000000 FALSE
15 0.8977333 0.0056862 -0.0022667 0.0000000 FALSE
20 0.8973000 0.0048474 -0.0027000 0.0000000 FALSE
30 0.8984667 0.0031528 -0.0015333 0.0000000 FALSE
50 0.9012400 0.0016474 0.0012400 0.0000001 FALSE
60 0.9033000 0.0013368 0.0033000 0.0000004 FALSE
100 0.8999400 0.0008495 -0.0000600 0.0000275 FALSE
200 0.9000500 0.0004848 0.0000500 0.0001946 FALSE
500 0.9007840 0.0002049 0.0007840 0.2189212 TRUE

Con la proporción de enfermos del 90% vemos que los histogramas al inicio no son muy simetricos, casi igual que el comportamiento de 10% con la diferencia que la prueba de shapiro no se acerco a considerar a ser normal, pero va subiendo el valor p.

Problema 4: Bootstrap

Se nos presenta un problema con una distribución desconocida (no parametrica), con la cual debemos aplicar una metodologia nueva de estimación de los parametros. el problema consta de el estudio de la medición de consumo de gasolina medido en millas/galon. debemos estimar un intervalo de confianza el consumo con un 95% confianza. para este ejercicio nos piden extraer \(k=1000\) de una pequeña muestra de datos:

1. Datos y estimados iniciales con bootstrap

Acontinuación vemos los datos suministrados

x=c( 7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)

Para hacer el bootstrap debemos tener los k eventos, y al ser 7 datos, podemos realizar una matriz de 1000 registros por estos 7 datos, para hacer la selección de los mismos empleariamos el metodo de reemplazo.(Es decir, el dato una vez sacado, puede volver a ser sacado)

boot=sample(x,1000*7,replace=TRUE) 
boot_matriz=matrix(boot,nrow=1000,ncol=7)

La estimador de cada muestra generada (o fila en este caso) sera la media

media_x=apply(boot_matriz,1,mean)

2. Calculando el intervalo de confianza IC

Se nos brindan dos metodos para calculoar mi intervalo de confianza con 95% confianza:

\[Método 1: (P_{2.5};P_{97.5})\] Que vendria siendo la identificación de los perceptiles 2.5 y y 97.5

ic_metodo1=quantile(media_x, probs=c(0.025, 0.975))
print(ic_metodo1)
##     2.5%    97.5% 
## 4.756929 6.460250

\[Método 2: (2\bar{X}−P_{97.5};2\bar{X}−P_{2.5})\] En este metodo empleamos un ajuste de los intervaos con los perceptiles 2.5 y 97.5

ic_metodo2=c(2*mean(media_x)-ic_metodo1[2], 2*mean(media_x)-ic_metodo1[1])
print(ic_metodo2)
##    97.5%     2.5% 
## 4.628293 6.331614

3. Grafico de los intervalos

hist(media_x, las=1, main=" ", ylab = " ", xlab = " ", col="#034A94")
abline(v=ic_metodo1, col="#FF7F00",lwd=2)
abline(v=ic_metodo2, col="#0EB0C6",lwd=2)

Se observa que le ejectivamente se puede esperar que el valor esperado este dentro de mi intervalo de confianza. De igual manera, se tiene que el metodo dos,(que tuvo una mejora a partir del metodo 1), tiene una mejor precisión y esta mas ajustado que su predecesor.