Problema 1

Introducción

A través del siguiente ejercicio, se estima el valor de \(\pi\) a partir de una simulación de ubicación de puntos en un plano espacial demarcado por un círculo de área \(\pi/4\).

En un primer momento, se hace la generación de n coordenadas valor mínimo de 0 y valor máximo de 1 con una distribución uniforme.

x <- runif(1000,0,1)

y <- runif(1000,0,1)

Así, teniendo en cuenta que resultan unas coordenadas \((Xi,Yi)\), se identifican los puntos que se ubican dentro del círculo a partir de la distancia que presentan en relación al centro del mismo, la cual debería ser menor a 0.5. De esta forma, se aplica la fórmula \((Xi-0.5)^2+(Yi-0.5)^2\)

distancia <- numeric() 
for (i in 1:length(x))
{ distancia[i] =
(x[i]-0.5)^2+(y[i]-0.5)^2
}

En este sentido, se asignan valores de 0 y 1 a los puntos que están fuera y dentro del círculo, respectivamente, para así definir la cantidad de puntos totales dentro del mismo:

puntos = ifelse(distancia < 0.25, 1, 0)
puntos_adentro<-sum(puntos)
print(paste0("Cantidad de puntos dentro del círculo: ", puntos_adentro))
## [1] "Cantidad de puntos dentro del círculo: 797"

Con este resultado, se procede a estimar el valor de \(\pi\) y el error.

#¿Cuál es su estimación de π?
pi_1<-4 * sum(puntos)/length(x)
print(paste0("Estimación de π: ", pi_1))
## [1] "Estimación de π: 3.188"
#estimacion del error con 1000 puntos 
error_1<-4 * sum(puntos)/length(x) - pi

print(paste0("Error: ", error_1))
## [1] "Error: 0.0464073464102071"

Una vez realizado este ejercicio, se procede a repetirlo con 10000 y 100000 puntos.

10000 puntos:

#estimacion con 10000 puntos

x2 <- runif(10000,0,1) 
y2 <- runif(10000,0,1)

#determinar las distancias con 10000
distancia_2 <- numeric ()

for (i in 1:length(x2)) 
{ distancia_2[i] =
(x2[i]-0.5)^2+(y2[i]-0.5)^2
}

puntos_2 = ifelse(distancia_2 < 0.25, 1, 0)

#¿Cuántos de los puntos están dentro del círculo con 10000?
puntos_a10000<-sum(puntos_2)
print(paste0("Cantidad de puntos dentro del círculo: ", puntos_a10000))
## [1] "Cantidad de puntos dentro del círculo: 7837"
#¿Cuál es su estimación de π con 10000?
pi_2<-4 * sum(puntos_2)/length(x2)
print(paste0("Estimación de π: ", pi_2))
## [1] "Estimación de π: 3.1348"
#estimacion del error con 10000 puntos 
error_2<-4 * sum(puntos_2)/length(x2) - pi
print(paste0("Error ", error_2))
## [1] "Error -0.00679265358979331"

100000 puntos:

#estimacion con 100000 puntos

x3 <- runif(100000,0,1) 
y3 <- runif(100000,0,1)

#determinar las distancias con 100000
distancia_3 <- numeric ()

for (i in 1:length(x3)) 
{ distancia_3[i] =
(x3[i]-0.5)^2+(y3[i]-0.5)^2
}

puntos_3 = ifelse(distancia_3 < 0.25, 1, 0)

#¿Cuántos de los puntos están dentro del círculo con 100000?
puntosa_100000<-sum(puntos_3)
print(paste0("Cantidad de puntos dentro del círculo: ", puntos_a10000))
## [1] "Cantidad de puntos dentro del círculo: 7837"
#¿Cuál es su estimación de π con 10000?
pi_3<-4 * sum(puntos_3)/length(x3)
print(paste0("Estimación de π: ", pi_3))
## [1] "Estimación de π: 3.13668"
#estimacion del error con 10000 puntos 
error_3<-4 * sum(puntos_3)/length(x3) - pi
print(paste0("Error ", error_3))
## [1] "Error -0.00491265358979298"

Como se puede ver, el error disminuye a medida que aumenta la cantidad de puntos en el estudio.

Problema 2

Introducción

El presente documento adelanta un ejercicio de simulación a partir de la generación de valores en una distribución exponencial. Su principal objetivo es la evaluación de cuatro (4) estimadores según su insesgadez, consistencia y eficiencia. En este sentido, los datos de punto de partida son:

Parametro θ: desconocido, se asigna en el ejercicio.

Tamaño de muestras: 4

Cantidad de experimentos: 20, 50, 100 y 1000

Estimadores a evaluar:

\[ a) \widehat{\theta_1}=\frac{(X_1+X_2)}{6}+\frac{(X_3+X_4)}{3} \]

\[b) \widehat{\theta_2}= \frac{(X_1+2X_2+3X_3+4X_4)}{5}\] \[ c) \widehat{\theta_3}= \frac {X_1+X_2+X_3+X_4}{4}\]

\[ d) \widehat{\theta_4}= \frac {min(X_1, X_2,X_3, X_4)+max(X_1, X_2, X_3, X_4)}{2}\]

#Creación de Funciones de estimadores
f1<-function(x) {((x[1]+x[2])/6)+((x[3]+x[4])/3)}
f2<-function(x) {(x[1]+(2*x[2])+(3*x[3])+(4*x[4]))/5}
f3<-function(x) {(x[1]+x[2]+x[3]+x[4])/4}
f4<-function(x) {(min(x)+max(x))/2}

# definición de tamaños de muestras
muestra=4
n1<-20 
n2<-50
n3<-100
n4<-1000
#definición de parametro θ
theta<-2
# librerias
library(ggplot2)

Evaluación con 20 iteraciones

#semilla
set.seed(123)
# m tamaño de replicas del experimento 
m1=n1*muestra  
# matriz de datos m x n 
sim1=matrix(rexp(m1, rate=theta), ncol=muestra)
# cálculo de los m valores para el estimador 1 
est1_20<-apply(sim1, 1,f1)
# cálculo de los m valores para el estimador 2
est2_20<-apply(sim1,1,f2) 
# cálculo de los m valores para el estimador 3
est3_20<-apply(sim1,1,f3) 
# cálculo de los m valores para el estimador 4
est4_20<-apply(sim1,1,f4)  
# data.frame con los tres estimadores data de m x 3
estimadores_20=data.frame(est1_20,est2_20,est3_20,est4_20) 

#Creación de matrices para análisis conjunto final
medias_20<- apply(estimadores_20,2, mean)
tmedias_20<- matrix(c(medias_20), ncol=4, byrow=TRUE)

sesgos_20 <- medias_20 - theta
tsesgos_20<- matrix(c(sesgos_20), ncol=4, byrow=TRUE)

eficiencias_20 <- apply(estimadores_20,2, var)
teficiencias_20<- matrix(c(eficiencias_20), ncol=4, byrow=TRUE)
#Generación de boxplot
label=c("Estimador 1","Estimador 2","Estimador 3", "Estimador 4")
boxplot(estimadores_20, main="Estimadores con n=20", names=label)
abline(h=theta,  col='red')

Como se puede ver, el estimador 2 presenta valores más cercanos al parámetro θ y los tres restantes con resultados similares relativamente más alejados del parámetro.

Evaluación con 50 iteraciones

#semilla
set.seed(123)
# m tamaño de replicas del experimento
m2=n2*muestra
# matriz de datos m x n
sim2=matrix(rexp(m2, rate=theta), ncol=muestra) 
# cálculo de los m valores para el estimador 1 
est1_50<-apply(sim2, 1,f1)
# cálculo de los m valores para el estimador 2
est2_50<-apply(sim2,1,f2) 
# cálculo de los m valores para el estimador 3
est3_50<-apply(sim2,1,f3) 
# cálculo de los m valores para el estimador 4
est4_50<-apply(sim2,1,f4)  
# data.frame con los tres estimadores data de m x 3
estimadores_50=data.frame(est1_50,est2_50,est3_50,est4_50) 

#Creación de matrices para análisis conjunto final
medias_50<- apply(estimadores_50,2, mean)
tmedias_50<- matrix(c(medias_50), ncol=4, byrow=TRUE)

sesgos_50 <- medias_50 - theta
tsesgos_50<- matrix(c(sesgos_50), ncol=4, byrow=TRUE)

eficiencias_50 <- apply(estimadores_50,2, var)
teficiencias_50<- matrix(c(eficiencias_50), ncol=4, byrow=TRUE)
#Generación de boxplot
label=c("Estimador 1","Estimador 2","Estimador 3", "Estimador 4")
boxplot(estimadores_50, main="Estimadores con n=50", names=label)
abline(h=theta,  col='red')

Al igual que las carácterísticas observadas en la realización de 20 experimentos, cuando se realizan 50 los resultados arrojan una mayor precisión del estimador 2 frente a su cercanía al valor del parámetro θ. En este sentido, dicho patrón se repite en la realización de 100 y 1000 experimentos:

#semilla
set.seed(123)
# m tamaño de replicas del experimento
m3=n3*muestra
# matriz de datos m x n
sim3=matrix(rexp(m3, rate=theta), ncol=muestra) 
# cálculo de los m valores para el estimador 1 
est1_100<-apply(sim3, 1,f1)
# cálculo de los m valores para el estimador 2
est2_100<-apply(sim3,1,f2) 
# cálculo de los m valores para el estimador 3
est3_100<-apply(sim3,1,f3) 
# cálculo de los m valores para el estimador 4
est4_100<-apply(sim3,1,f4)  
# data.frame con los tres estimadores data de m x 3
estimadores_100=data.frame(est1_100,est2_100,est3_100,est4_100) 

#Creación de matrices para análisis conjunto final
medias_100<- apply(estimadores_100,2, mean)
tmedias_100<- matrix(c(medias_100), ncol=4, byrow=TRUE)

sesgos_100 <- medias_100 - theta
tsesgos_100<- matrix(c(sesgos_100), ncol=4, byrow=TRUE)

eficiencias_100 <- apply(estimadores_100,2, var)
teficiencias_100<- matrix(c(eficiencias_100), ncol=4, byrow=TRUE)
#Generación de boxplot
label=c("Estimador 1","Estimador 2","Estimador 3", "Estimador 4")
boxplot(estimadores_100, main="Estimadores con n=100", names=label)
abline(h=theta,  col='red') 

#semilla
set.seed(123)
# m tamaño de replicas del experimento
m4=n4*muestra
# matriz de datos m x n
sim4=matrix(rexp(m4, rate=theta), ncol=muestra) 
# cálculo de los m valores para el estimador 1 
est1_1000<-apply(sim4, 1,f1)
# cálculo de los m valores para el estimador 2
est2_1000<-apply(sim4,1,f2) 
# cálculo de los m valores para el estimador 3
est3_1000<-apply(sim4,1,f3) 
# cálculo de los m valores para el estimador 4
est4_1000<-apply(sim4,1,f4)  
# data.frame con los tres estimadores data de m x 3
estimadores_1000=data.frame(est1_1000,est2_1000,est3_1000,est4_1000) 

#Creación de matrices para análisis conjunto final
medias_1000<- apply(estimadores_1000,2, mean)
tmedias_1000<- matrix(c(medias_1000), ncol=4, byrow=TRUE)

sesgos_1000 <- medias_1000 - theta
tsesgos_1000<- matrix(c(sesgos_1000), ncol=4, byrow=TRUE)

eficiencias_1000 <- apply(estimadores_1000,2, var)
teficiencias_1000<- matrix(c(eficiencias_1000), ncol=4, byrow=TRUE)
#Generación de boxplot
label=c("Estimador 1","Estimador 2","Estimador 3", "Estimador 4")
boxplot(estimadores_1000, main="Estimadores con n=1000", names=label)
abline(h=theta,  col='red')

Evaluación de características

Medias (para analizar consistencia)

tab_medias<-rbind(medias_20, medias_50, medias_100, medias_1000)
colnames(tab_medias) <- c('Estimador 1', 'Estimador 2', 'Estimador 3', 'Estimador 4')
rownames(tab_medias) <- c('Medias - 20', 'Medias - 50', 'Medias - 100',
                   'Medias - 1000')
print(tab_medias)
##               Estimador 1 Estimador 2 Estimador 3 Estimador 4
## Medias - 20     0.5177122   1.0219794   0.5085644   0.6003939
## Medias - 50     0.4972023   0.9895097   0.5036166   0.5788121
## Medias - 100    0.4941755   0.9783242   0.4965358   0.5733920
## Medias - 1000   0.4980205   0.9898509   0.4999876   0.5810203

Como se puede observar, en general, según lo señalado anteriormente, el estimador 2 es el que más tiende a acercarse al parámetro θ. De manera que, si bien la media de todos los estimadores tiende a aumentar para así acercarse al valor de theta (2), sólo el estimador 2 muestra una real aproximación a dicho valor.

Sesgos

tab_sesgos<-rbind(sesgos_20, sesgos_50,
           sesgos_100, sesgos_1000)
colnames(tab_sesgos) <- c('Estimador 1', 'Estimador 2', 'Estimador 3', 'Estimador 4')
rownames(tab_sesgos) <- c('Sesgos - 20', 'Sesgos - 50', 'Sesgos - 100','Sesgos - 1000')
tab_sesgos<-as.data.frame(tab_sesgos)
print(tab_sesgos)
##               Estimador 1 Estimador 2 Estimador 3 Estimador 4
## Sesgos - 20     -1.482288  -0.9780206   -1.491436   -1.399606
## Sesgos - 50     -1.502798  -1.0104903   -1.496383   -1.421188
## Sesgos - 100    -1.505824  -1.0216758   -1.503464   -1.426608
## Sesgos - 1000   -1.501980  -1.0101491   -1.500012   -1.418980

En el caso de la insesgadez, el 2 estimador presenta mejores características, ya que para las 4 cantidades de experimentos planteados, el valor esperado de dicho estimador presenta menores diferencias frente al parámetro θ.

Eficiencias

tab_efic<-rbind(eficiencias_20, eficiencias_50,
                eficiencias_100, eficiencias_1000)
colnames(tab_efic) <- c('Estimador 1', 'Estimador 2', 'Estimador 3', 'Estimador 4')
rownames(tab_efic) <- c('Varianza - 20', 'Varianza - 50', 'Varianza - 100','Varianza - 1000')
print(tab_efic)
##                 Estimador 1 Estimador 2 Estimador 3 Estimador 4
## Varianza - 20    0.06803333   0.2192192  0.04964954   0.1301817
## Varianza - 50    0.06697232   0.2679561  0.06724301   0.1232753
## Varianza - 100   0.05180615   0.2348037  0.04945630   0.0822846
## Varianza - 1000  0.07155260   0.2941633  0.06496121   0.1009200

No obstante, al estudiar las eficiencias se observa que, a pesar de presentar valores más cercanos a theta, el estimador 2 es el de mayores varianzas en los 4 casos, siendo el estimador 3 el más eficiente seguido del 1.

Conclusión

Si bien, en principio, el estimador 2 presenta valores cercanos al parámetro a lo largo de k=20, 50, 100 y 1000, no revela la misma eficiencia de los otros estimadores estudiados, razón por la cual se presentan dificultades en la definición de un mejor estimador en función del parámetro del modelo.

Problema 3

Introducción

El siguiente ejercicio, plantea la ejecución de pasos específicos para la verificación del Teorema de Límite Central:

a) Simulación para generar una población de n=1000 (Lote), donde el porcentaje de individuos (supongamos plantas) enfermas sea del 50%.

G_valores <- function (n,p_sanas,p_enfermas){
estado <- c(rep("sanas",n*p_sanas), rep("enfermas",n*p_enfermas))
print(paste("lote:",n,"sanas:", n*p_sanas,"enfermas:", n*p_enfermas))
return(estado)
}
n=1000
p_sanas=0.5
p_enfermas=0.5
lote<- G_valores(n,p_sanas,p_enfermas)
## [1] "lote: 1000 sanas: 500 enfermas: 500"

b) Generación de una función que permite obtener una muestra aleatoria de la población y cálculo del estimador de la proporción muestral \(\hat{p}\) para un tamaño de muestra dado \(n\).

generar_Muestra <- function(lote, n) {
  muestra <- sample(lote, n)
  return(muestra)
}
calcularProporcion <- function(lote, n, estado) {
  muestra <- generar_Muestra(lote, n)
  P <- sum(muestra == estado) / n
  return(P)
}
proporcion <- calcularProporcion(lote, 100, "enfermas")

print(paste("Para una muestra de tamaño: ", 100, "se obtuvo un ^P =",proporcion))
## [1] "Para una muestra de tamaño:  100 se obtuvo un ^P = 0.62"

c) En este sentido, se repite el escenario anterior 500 veces y se analizan los resultados en cuanto al comportamiento de los 500 resultados del estimador \(\hat{p}\). 

tam_muestra<- 100 
repeticiones<-500
estimaciones_proporcion <- numeric(repeticiones) 

for (i in 1:repeticiones) {
  estimaciones_proporcion[i] <-calcularProporcion(lote, tam_muestra,"enfermas")
} # Se realizo muestreo  y se registro las estimaciones
hist(estimaciones_proporcion,col="#377eb8", main ="Histograma 500 repeticiones",
     xlab="Probabilidad de las muestras", ylab="Frecuencia", las=1,
     font.axis=4)
line = mean(estimaciones_proporcion)
abline (v=line, lwd = 4, lty = 2, col="darkred")

media <- mean.default(estimaciones_proporcion, na.rm = TRUE)
mediana <- median.default(estimaciones_proporcion, na.rm = TRUE)
var <- var(estimaciones_proporcion, na.rm = TRUE)
desvest <- sd(estimaciones_proporcion, na.rm = TRUE)
q1 <- quantile(estimaciones_proporcion, probs = 0.25, na.rm = TRUE)
q3 <- quantile(estimaciones_proporcion, probs = 0.75, na.rm = TRUE)
print(paste("Media:", media))
## [1] "Media: 0.49982"
print(paste("Mediana:", mediana))
## [1] "Mediana: 0.5"
print(paste("Varianza:", var))
## [1] "Varianza: 0.00227511783567134"
print(paste("Desviación Estándar:", desvest))
## [1] "Desviación Estándar: 0.0476981953083274"
print(paste("Cuartil 1:", q1))
## [1] "Cuartil 1: 0.47"
print(paste("Cuartil 3:", q3))
## [1] "Cuartil 3: 0.53"

¿Qué tan simétricos o sesgados son los resultados obtenidos? y ¿qué se puede observar en cuanto a la variabilidad?

En el historgrama se observa que, tomando una muestra de 100 datos, se tiene que el conjunto de estimadores obtenidos bajo las 500 iteraciones presentan un comportamiento de una distribución normal ya que el grafico muestra una similitud con la campana de Gauss, exponiendo alta simetría. Como se observa en el gráfico anterior, la media (0.497) y la mediana (0.495) son muy similares y concuerdan con la proporción poblacional. Ademas, revisando el calculo de la varianza, que es una medida de dispersion (0.0022134), se observa que tiene un valor tendiente a cero, por ende, se trata de una muestra de datos confiable donde no hay datos que la sesguen.

d) Se repiten los puntos b y c para tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500 y se comparan los resultados obtenidos para los diferentes tamaños de muestra en cuanto a la normalidad. Se aplican pruebas de bondad y ajuste (shapiro wilks :shspiro.test()) y métodos gráficos (gráfico de normalidad: qqnorm()). 

item <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
for(i in item) {
  resultados <- replicate(n, calcularProporcion(lote,i, "enfermas"))
  par(mfrow=c(1,2))
  hist(resultados, col="#377eb8", main = paste("Histograma de la muestra n =", i),     xlab="Probabilidad de la muestra", ylab="Frecuencia", las=1, font.axis=4)
  abline (v=mean(resultados), lwd = 4, lty = 2, col="#ff7f00")
  shapiroRes <- shapiro.test(resultados)
  print(shapiroRes)
  qqnorm(resultados, col = "#377eb8")
  qqline(resultados, col = "#ff7f00")

}
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.92929, p-value < 2.2e-16

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.96267, p-value = 2.59e-15

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97527, p-value = 5.117e-12

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97847, p-value = 5.261e-11

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.98668, p-value = 6.774e-08

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.99029, p-value = 3.585e-06

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.99217, p-value = 3.817e-05

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.99523, p-value = 0.00321

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.99686, p-value = 0.04529

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.99692, p-value = 0.05035

En este caso se evidencia que no todas las muestras siguen una distribucion normal, ya que cuando n = 5,10,15,30,50 su valor \(p\) toma un valor muy pequeno, de manera que se rechaza la hipotesis nula de normalidad. En el caso de las muestras n>50 el valor \(p\) es ligeramente superior a 0.05, lo que indica que no hay suficiente evidencia para rechazar la hipótesis nula de normalidad a un nivel de significancia del 5%, considerándose como un margen estrecho. Adicional, se observa en los histogramas los comportamientos asimétricos en las simulaciones iniciales que convergen hacia la campana de Gauss a medida que incrementa \(n\).

Por otra parte, revisando el comportamiento de sus correspondientes gráficos Q-Q, se observa que la línea de regresión de las coordenadas entre cuantiles téoricos y empíricos tiende a acercarse a medida que aumenta la muestra de los datos.

e) Se repite toda la simulación (puntos a – d), pero ahora para lotes con 10% de plantas enfermas y de nuevo para lotes con un 90% de plantas enfermas.

Lotes con 10% de plantas enfermas

lote <- G_valores(n=1000,p_sanas=0.9,p_enfermas=0.1)
## [1] "lote: 1000 sanas: 900 enfermas: 100"
item <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
for(n in item) {
  resultados <- replicate(n, calcularProporcion(lote, n, "enfermas"))
  par(mfrow=c(1,2))
  
hist(resultados, col="#377eb8", main = paste("Histograma de la muestra n =", n), xlab="Probabilidad de las muestras", ylab="Frecuencia", las=1, font.axis=4)
  abline (v=mean(resultados), lwd = 4, lty = 2, col="#ff7f00")
  shapiroRes <- shapiro.test(resultados)
  print(shapiroRes)
  qqnorm(resultados, col = "#377eb8")
  qqline(resultados, col = "#ff7f00")
  
}
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.68403, p-value = 0.00647

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.79406, p-value = 0.01228

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.88143, p-value = 0.04983

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.92581, p-value = 0.1282

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.90033, p-value = 0.008559

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.96311, p-value = 0.12

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97545, p-value = 0.2669

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.98208, p-value = 0.1924

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.98407, p-value = 0.02321

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.99451, p-value = 0.07015

Lotes con un 90% de plantas enfermas

lote <- G_valores(n=1000,p_sanas=0.1,p_enfermas=0.9)
## [1] "lote: 1000 sanas: 100 enfermas: 900"
item <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
for(n in item) {
  resultados <- replicate(n, calcularProporcion(lote, n, "enfermas"))
  par(mfrow=c(1,2))
  
hist(resultados, col="#377eb8", main = paste("Histograma de la muestra n =", n), xlab="Probabilidad de las muestras", ylab="Frecuencia", las=1, font.axis=4)
  abline (v=mean(resultados), lwd = 4, lty = 2, col="#ff7f00")
  shapiroRes <- shapiro.test(resultados)
  print(shapiroRes)
  qqnorm(resultados, col = "#377eb8")
  qqline(resultados, col = "#ff7f00")
  
}
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.68403, p-value = 0.00647

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.90444, p-value = 0.2449

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.8959, p-value = 0.0824

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.83373, p-value = 0.00288

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.90751, p-value = 0.01288

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.93594, p-value = 0.009336

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97332, p-value = 0.2114

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.98372, p-value = 0.2559

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.98867, p-value = 0.1134

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.99279, p-value = 0.0166

Realizando la comparacion de los dos casos se observa que no en todas las muestras los datos siguen una distribucion normal; en el caso en donde el 90% de las plantas están enfermas se rechaza la hipotesis nula de distribucion normal para las muestras: 5, 10, 20, 100 y 200, ya que su valor \(p\) es menor a 0.05 y significa que la hipótesis nula es falsa. Para el caso en donde el 10% de las plantas está enferma las muestras: 10, 20,30, 60,200 y 500 hay un p > 0,05 que indica que la hipótesis nula es verdadera y, por ende, sigue una distribucion normal. Esto tambien se observa en los gráficos, ya que proporcionaron una representación visual para evaluar la normalidad de los datos.

Por ultimo, además, se observó un patró de sesgo direccional de los datos,en donde una proporción del 10% de plantas enferma resultó en sesgo hacia la derecha y una proporción del 90% de plantas enferma generó sesgo hacia la izquierda, lo que refleja la concentración de datos alrededor de la proporción del parámetro.

Problema 4

Introducción

A continuación, se realiza un ejercicio de simulación de muestras Bootstrap a partir de los siguientes datos, que corresponden a la medición de la eficiencia de combustible (millas/galón) en 7 camiones: 7.69, 4.97, 4.56, 6.49, 4.34, 6.24 y 4.45.

A partir de los mismos, se realizan 1000 experimentos y se busca construir los siguientes intervalos de confianza conformados por los percentiles P2.5 y P97.5:

Método 1: \[ P_{2.5}; P_{97.5}\]

Método 2:

\[ (2\overline{X}−P_{97.5};2\overline{X}−P_{2.5}) \]

#Creación de matriz de muestras
# datos muestra
x=c( 7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45) 
# semilla reproducibilidad
set.seed(123)  
# se extraen n x m muestras
boot=sample(x,7000,replace=TRUE)   
# se construye matriz de n x m 
matriz_muestras=matrix(boot,ncol=7, byrow=TRUE)  
# se calculan m medias por fila
medias_muestras=apply(matriz_muestras,1,mean)  
medias_muestrasd=as.data.frame(medias_muestras)
print(head(matriz_muestras))
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 4.45 4.45 4.56 6.24 4.56 4.97 4.97
## [2,] 6.24 4.56 4.34 6.49 6.24 6.24 7.69
## [3,] 4.97 4.56 4.34 4.56 4.56 7.69 6.49
## [4,] 7.69 7.69 4.34 4.56 4.97 4.45 4.97
## [5,] 7.69 6.24 4.56 6.49 6.24 7.69 4.56
## [6,] 4.45 4.34 6.49 4.45 4.97 4.34 4.45

Así, se muestra la estructura de los datos resultantes del ejercicio de simulación y frente a los cuales se definieron las medias para el análisis de los cuartiles.

#Definición de cuartiles
# se calcula IC método 1
ic1=quantile(medias_muestras, probs=c(0.025, 0.975)) 
# se calcula IC método 2
ic2=c(2*mean(medias_muestras)-ic1[2], 2*mean(medias_muestras)-ic1[1]) 

tabla_cuant<- rbind(ic1, ic2)
rownames(tabla_cuant) =c('Método 1','Método 2')
print(tabla_cuant)
##              2.5%    97.5%
## Método 1 4.748393 6.508643
## Método 2 4.549526 6.309776

Cómo se puede ver, ambos métodos arrojan valores muy similares para los dos cuartiles, lo cual implica que tienen resultados, prácticamente, idénticos:

#Análisis de Intervalos de Confianza
valor_esperado<-apply(medias_muestrasd,2,FUN=mean)
media_datos=mean(x)
comparacion_medias<- rbind(media_datos, valor_esperado)
rownames(comparacion_medias) =c('Datos originales','Valor esperado')
colnames(comparacion_medias) = 'Media'

Además, resulta interesante resaltar que la media de los datos originales es bastante cercana al valor esperado resultante de los 1000 experimentos. Lo cual, gráficamente, se puede evidencia de la siguiente manera:

hist(medias_muestras, las=1, main="Medias de las Muestras Bootstrap", ylab = "Frecuencia", xlab = "Medias", col="#43766C", border="grey")
abline(v=ic1, col="orange", lwd=2)  # rojo para método 1
abline(v=ic2, col="purple", lwd=2)  # morado para método 2
abline(v=mean(x), col="red", lwd=3)
abline(v=valor_esperado, col="lightblue", lwd=3)
legend(6.3,195, # posición de la leyenda
       legend=c("IC Método 1", "IC Método 2", "Media de valores", "Valor Esperado"), # texto de las etiquetas
       col=c("orange", "purple","red","lightblue"),
       cex=0.7, # colores de las líneas
       lwd=2) # grosor de las líneas

Finalmente, se observa que la cantidad de resultados que se delimitan en los dos cuartiles, para ambos métodos es exactamente la misma, siendo 950 dicha cantidad (25 valores por fuera de cada límite para ambos casos), cumpliendo con el 95% del intervalo de confianza definido en el ejericio, haciéndolo confiable.

#cantidad de valores dentro de los intervalos
d1=(sum(medias_muestras>ic1[1]&medias_muestras<ic1[2])/1000) #método 1
d2=(sum(medias_muestras>ic2[1]&medias_muestras<ic2[2])/1000) #método 2