PUNTO 2 TALLER 1

#LIBRERIAS
library(plotly)

## Warning: package 'plotly' was built under R version 4.2.3

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 4.2.3

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(ggplot2)
library(ape)

## Warning: package 'ape' was built under R version 4.2.3

library(gstat)

## Warning: package 'gstat' was built under R version 4.2.3

library(distrEx)

## Warning: package 'distrEx' was built under R version 4.2.3

## Loading required package: distr

## Warning: package 'distr' was built under R version 4.2.3

## Loading required package: startupmsg

## Warning: package 'startupmsg' was built under R version 4.2.3

## Utilities for Start-Up Messages (version 0.9.6.1)

## For more information see ?"startupmsg", NEWS("startupmsg")

## Loading required package: sfsmisc

## Object Oriented Implementation of Distributions (version 2.9.3)

## Attention: Arithmetics on distribution objects are understood as operations on corresponding random variables (r.v.s); see distrARITH().
## Some functions from package 'stats' are intentionally masked ---see distrMASK().
## Note that global options are controlled by distroptions() ---c.f. ?"distroptions".

## For more information see ?"distr", NEWS("distr"), as well as
##   http://distr.r-forge.r-project.org/
## Package "distrDoc" provides a vignette to this package as well as to several extension packages; try vignette("distr").

## 
## Attaching package: 'distr'

## The following objects are masked from 'package:stats':
## 
##     df, qqplot, sd

## Extensions of Package 'distr' (version 2.9.2)

## Note: Packages "e1071", "moments", "fBasics" should be attached /before/ package "distrEx". See distrExMASK().Note: Extreme value distribution functionality has been moved to
##       package "RobExtremes". See distrExMOVED().

## For more information see ?"distrEx", NEWS("distrEx"), as well as
##   http://distr.r-forge.r-project.org/
## Package "distrDoc" provides a vignette to this package as well as to several related packages; try vignette("distr").

## 
## Attaching package: 'distrEx'

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, median, var

library(gridExtra)

## Warning: package 'gridExtra' was built under R version 4.2.3

cc = 1000376863
ultimo_digito <- as.numeric(substr(cc, nchar(cc), nchar(cc)))
if (!is.na(ultimo_digito)) {
  if (ultimo_digito == 1 || ultimo_digito == 2) {
    a <- 6
    b <- 24
  } else if (ultimo_digito == 3 || ultimo_digito == 4) {
    a <- 8
    b <- 18
  } else if (ultimo_digito == 5 || ultimo_digito == 4) {
    a <- 12
    b <- 12
  } else if (ultimo_digito == 6 || ultimo_digito == 7) {
    a <- 16
    b <- 9
  } else if (ultimo_digito == 8 || ultimo_digito == 9) {
    a <- 18
    b <- 8
  } else {
    a <- NA
    b <- NA
  }
} else {
  a <- NA
  b <- NA
}
mensaje <- paste("El tamaño de la grilla es", a, "x", b)
print(mensaje)

## [1] "El tamaño de la grilla es 8 x 18"

Creación de Datos y Grilla

resp = runif(n = 144, min = 0.66, max = 0.86)
xy = expand.grid(x = 1:a, y = 1:b)
plot(xy,cex=0.5,pch=15)

# UNION DE LOS DATOS EN DATA FRAME
datos1 = data.frame(resp, x = xy$x, y = xy$y)

Visualización de los datos

ggplot(datos1, aes(x = x, y = y, fill = resp)) +
  geom_tile() +
  theme_minimal() +
  scale_fill_gradient(low = "green", high = "red")

m_dist=as.matrix(dist(cbind(datos1$x,datos1$y)))
m_resp=1/m_dist
diag(m_resp)<-0
w=m_resp

Indice de MORAN 1

Moran.I(datos1$resp, w)

## $observed
## [1] -0.01472526
## 
## $expected
## [1] -0.006993007
## 
## $sd
## [1] 0.008379295
## 
## $p.value
## [1] 0.3561216

Creación de Simulaciones (10 RANGOS de 100 simulaciones)

NDVI_all = c() 
for(i in seq(0.1,5,0.5)){
  g.dummy = gstat(formula = z~1, locations = ~x+y,
                  dummy = T, beta = 0.5, model = vgm(psill = 0.0009,
                                                     range = i, model = 'Sph'), nmax = 20)
  NDVI = predict(g.dummy, newdata = xy, nsim = 100)
  NDVI_all = c(NDVI_all, NDVI)
}

## [using unconditional Gaussian simulation]
## [using unconditional Gaussian simulation]
## [using unconditional Gaussian simulation]
## [using unconditional Gaussian simulation]
## [using unconditional Gaussian simulation]
## [using unconditional Gaussian simulation]
## [using unconditional Gaussian simulation]
## [using unconditional Gaussian simulation]
## [using unconditional Gaussian simulation]
## [using unconditional Gaussian simulation]

datos = as.data.frame(NDVI_all)

Función para separar las tablas por cada Rango

tablas <- list()
n_columnas <- 102

for (i in 1:10) {
  start_col <- (i - 1) * n_columnas + 1
  end_col <- i * n_columnas
  sub_tabla <- datos[, start_col:end_col]
  #calculo de la media
  suma_columnas <- rowSums(sub_tabla[, 3:102], na.rm = TRUE)
  media_columnas <- suma_columnas / 100
  #añadir columna con la media al final de la tabla
  sub_tabla <- cbind(sub_tabla, media = media_columnas)
  tablas[[i]] <- sub_tabla
}

Creación de las tablas

r1 = tablas[[1]]
r2 = tablas[[2]]
r3 = tablas[[3]]
r4 = tablas[[4]]
r5 = tablas[[5]]
r6 = tablas[[6]]
r7 = tablas[[7]]
r8 = tablas[[8]]
r9 = tablas[[9]]
r10 = tablas[[10]]

Indice de Moran para cada rango

m1 <- Moran.I(r1$media, w)$observed
m2 <- Moran.I(r2$media, w)$observed
m3 <- Moran.I(r3$media, w)$observed
m4 <- Moran.I(r4$media, w)$observed
m5 <- Moran.I(r5$media, w)$observed
m6 <- Moran.I(r6$media, w)$observed
m7 <- Moran.I(r7$media, w)$observed
m8 <- Moran.I(r8$media, w)$observed
m9 <- Moran.I(r9$media, w)$observed
m10 <- Moran.I(r10$media, w)$observed
tablamoran = as.data.frame(rbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10))

Visualización de los Indices de Moran

barplot(tablamoran$V1)

# Creación de tratamientos
tratamientos = sample(rep(paste0("T", 1:12), each = 12))
print(tratamientos)

##   [1] "T10" "T11" "T5"  "T1"  "T7"  "T10" "T12" "T2"  "T3"  "T4"  "T8"  "T8" 
##  [13] "T9"  "T8"  "T12" "T3"  "T2"  "T12" "T3"  "T11" "T2"  "T5"  "T3"  "T9" 
##  [25] "T11" "T7"  "T4"  "T5"  "T3"  "T11" "T11" "T9"  "T5"  "T5"  "T12" "T12"
##  [37] "T1"  "T12" "T8"  "T6"  "T11" "T1"  "T7"  "T6"  "T9"  "T4"  "T10" "T6" 
##  [49] "T2"  "T9"  "T12" "T4"  "T10" "T4"  "T1"  "T2"  "T6"  "T6"  "T9"  "T11"
##  [61] "T2"  "T2"  "T12" "T7"  "T12" "T1"  "T4"  "T10" "T6"  "T9"  "T7"  "T6" 
##  [73] "T5"  "T8"  "T10" "T12" "T10" "T4"  "T12" "T8"  "T3"  "T10" "T5"  "T7" 
##  [85] "T9"  "T7"  "T7"  "T10" "T4"  "T9"  "T10" "T6"  "T11" "T5"  "T4"  "T6" 
##  [97] "T11" "T3"  "T7"  "T1"  "T8"  "T4"  "T2"  "T6"  "T5"  "T5"  "T3"  "T10"
## [109] "T2"  "T9"  "T7"  "T10" "T3"  "T8"  "T2"  "T4"  "T5"  "T1"  "T11" "T1" 
## [121] "T4"  "T11" "T11" "T9"  "T1"  "T12" "T8"  "T1"  "T2"  "T8"  "T7"  "T8" 
## [133] "T6"  "T1"  "T3"  "T3"  "T7"  "T2"  "T3"  "T8"  "T6"  "T9"  "T5"  "T1"

# Añadir tratamientos a tabla de rango
r1 <- cbind(Tratamientos = tratamientos, r1)

Visualización de la distribución de los tratamientos

ggplot(r1,aes(x=x,y=y,fill=sim1,label=tratamientos))+
  geom_tile()+
  geom_label(fill='white')+
  theme_minimal()+
theme(legend.position = "none")

# crear una submuestra con columnas de interes
sbt1 <- r1[, 4:103]
#ANOVA para cada columna con respecto a los tratamientos
ranova1 <- lapply(sbt1, function(x) {
aov(x ~ r1[[1]])})
#Extraer el estadístico F de cada análisis de varianza
ef1 <- sapply(ranova1, function(x) {
summary(x)[[1]][1, "F value"]})
# Mostrar los estadísticos F
print(ef1)

##      sim1      sim2      sim3      sim4      sim5      sim6      sim7      sim8 
## 2.1388343 2.2944939 1.7490072 1.0834997 1.5940397 0.4727554 0.9183339 0.7686974 
##      sim9     sim10     sim11     sim12     sim13     sim14     sim15     sim16 
## 1.2601671 1.1552079 0.5369350 0.2354232 0.7985152 2.4570446 1.1372988 0.5202359 
##     sim17     sim18     sim19     sim20     sim21     sim22     sim23     sim24 
## 0.6189676 0.8853988 1.2049997 1.0166760 0.8774161 0.9481890 0.7962218 1.5618160 
##     sim25     sim26     sim27     sim28     sim29     sim30     sim31     sim32 
## 0.7224322 1.1001901 1.2240019 1.5745246 1.0946127 0.7063647 0.6987952 1.2048476 
##     sim33     sim34     sim35     sim36     sim37     sim38     sim39     sim40 
## 0.9048847 0.5470431 1.1205362 0.3367344 0.7570170 0.7561494 0.7708435 0.8002546 
##     sim41     sim42     sim43     sim44     sim45     sim46     sim47     sim48 
## 2.2545901 1.1456483 0.6137046 0.7914772 1.1009087 0.7705305 0.4708138 0.9102053 
##     sim49     sim50     sim51     sim52     sim53     sim54     sim55     sim56 
## 1.2116078 0.4561146 0.7881033 0.3385155 0.5779626 2.0376076 1.0039071 0.8048103 
##     sim57     sim58     sim59     sim60     sim61     sim62     sim63     sim64 
## 0.7106881 1.1580947 0.5996784 1.1097986 0.8390127 1.6079663 0.8292246 0.4017533 
##     sim65     sim66     sim67     sim68     sim69     sim70     sim71     sim72 
## 1.2547950 0.6171318 1.0121941 1.1856063 1.6391598 1.0367400 1.4290615 0.7317164 
##     sim73     sim74     sim75     sim76     sim77     sim78     sim79     sim80 
## 1.4140651 1.2631775 1.4887665 0.8388356 0.5378556 1.0008751 0.8385581 0.9567163 
##     sim81     sim82     sim83     sim84     sim85     sim86     sim87     sim88 
## 0.4658013 1.5376430 1.2682707 1.9181415 0.5394461 0.8096610 0.8115995 0.2937064 
##     sim89     sim90     sim91     sim92     sim93     sim94     sim95     sim96 
## 1.9818996 1.1693243 1.1188562 0.9835255 1.4097381 0.8908283 1.4277795 1.1208333 
##     sim97     sim98     sim99    sim100 
## 0.9001482 2.0597613 1.0825013 0.6159270

# DATOS UNIVERSALES (GRADOS DE LIBERTAD)
gl_numerador = 11
gl_denominador = 132
cuantil_95 <- qf(0.95, gl_numerador, gl_denominador)

Graficas del F valor por rango

hist(ef1,
     main = "Estadísticos F Rango1",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef1), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

r2 <- cbind(Tratamientos = tratamientos, r2)

sbt2 <- r2[, 4:103]

  ranova2 <- lapply(sbt2, function(x) {
    aov(x ~ r2[[1]])})
  
  ef2 <- sapply(ranova2, function(x) {
    summary(x)[[1]][1, "F value"]})
  
  
hist(ef2,
     main = "Estadísticos F Rango2",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef2), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

r3 <- cbind(Tratamientos = tratamientos, r3)

sbt3 <- r3[, 4:103]

  ranova3 <- lapply(sbt3, function(x) {
    aov(x ~ r3[[1]])})
  
  ef3 <- sapply(ranova3, function(x) {
    summary(x)[[1]][1, "F value"]})

  
hist(ef3,
     main = "Estadísticos F Rango3",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef3), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

r4 <- cbind(Tratamientos = tratamientos, r4)

sbt4 <- r4[, 4:103]

  ranova4 <- lapply(sbt4, function(x) {
    aov(x ~ r4[[1]])})
  
  ef4 <- sapply(ranova4, function(x) {
    summary(x)[[1]][1, "F value"]})
  
  
hist(ef4,
     main = "Estadísticos F Rango4",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef4), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

r5 <- cbind(Tratamientos = tratamientos, r5)

sbt5 <- r5[, 4:103]

  ranova5 <- lapply(sbt5, function(x) {
    aov(x ~ r5[[1]])})
  
  ef5 <- sapply(ranova5, function(x) {
    summary(x)[[1]][1, "F value"]})
  
hist(ef5,
     main = "Estadísticos F Rango5",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef5), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

r6 <- cbind(Tratamientos = tratamientos, r6)

sbt6 <- r6[, 4:103]

  ranova6 <- lapply(sbt6, function(x) {
    aov(x ~ r6[[1]])})
  
  ef6 <- sapply(ranova6, function(x) {
    summary(x)[[1]][1, "F value"]})
  
hist(ef6,
     main = "Estadísticos F Rango6",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef6), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

r7 <- cbind(Tratamientos = tratamientos, r7)

sbt7 <- r7[, 4:103]

  ranova7 <- lapply(sbt7, function(x) {
    aov(x ~ r7[[1]])})
  
  ef7 <- sapply(ranova7, function(x) {
    summary(x)[[1]][1, "F value"]})
  
hist(ef7,
     main = "Estadísticos F Rango7",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef7), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

r8 <- cbind(Tratamientos = tratamientos, r8)

sbt8 <- r8[, 4:103]

  ranova8 <- lapply(sbt8, function(x) {
    aov(x ~ r8[[1]])})
  
  ef8 <- sapply(ranova8, function(x) {
    summary(x)[[1]][1, "F value"]})

hist(ef8,
     main = "Estadísticos F Rango8",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef8), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

r9 <- cbind(Tratamientos = tratamientos, r9)

sbt9 <- r9[, 4:103]

  ranova9 <- lapply(sbt9, function(x) {
    aov(x ~ r9[[1]])})
  
  ef9 <- sapply(ranova9, function(x) {
    summary(x)[[1]][1, "F value"]})
  
hist(ef9,
     main = "Estadísticos F Rango9",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef9), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

r10 <- cbind(Tratamientos = tratamientos, r10)

sbt10 <- r10[, 4:103]

  ranova10 <- lapply(sbt10, function(x) {
    aov(x ~ r10[[1]])})
  
  ef10 <- sapply(ranova10, function(x) {
    summary(x)[[1]][1, "F value"]})
  
hist(ef10,
     main = "Estadísticos F Rango10",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)
abline(v = mean(ef10), col = "red", lwd = 2)
abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE)

Graficar F valores

# Configurar la disposición de la ventana gráfica para 2 filas y 5 columnas
par(mfrow = c(2, 5))

# Suponiendo que tienes 10 conjuntos de datos ef1, ef2, ..., ef10 y cuantil_95_1, cuantil_95_2, ..., cuantil_95_10

for (i in 1:10) {
  # Genera los nombres de las variables dinámicamente
  ef <- get(paste0("ef", i))
  
  # Crear el histograma
  hist(ef,
       main = paste("Rango", i),
       xlab = "Valores del Estadístico F",
       ylab = "Frecuencia",
       col = "blue",
       border = "black",
       breaks = 50)
  
  # Añadir líneas verticales y curva de distribución
  abline(v = mean(ef), col = "red", lwd = 2)
  abline(v = cuantil_95, col = "red", lwd = 2, lty = 2)
  curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
        col = "darkgreen", 
        lwd = 2, 
        add = TRUE)
}

# Restablecer la disposición de la ventana gráfica
par(mfrow = c(1, 1))

Graficas solo cuantil 95 por rango

# Filtrar los valores mayores o iguales al cuantil 95
ef1_95 <- ef1[ef1 >= cuantil_95]

# Crear el histograma solo con los valores filtrados
hist(ef1_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango1",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

# Añadir la curva de la distribución F solo para los valores mayores o iguales al cuantil 95
curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

ef2_95 <- ef2[ef2 >= cuantil_95]

hist(ef2_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango2",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

ef3_95 <- ef3[ef3 >= cuantil_95]

hist(ef3_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango3",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

ef4_95 <- ef4[ef4 >= cuantil_95]

hist(ef4_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango4",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

ef5_95 <- ef5[ef5 >= cuantil_95]

hist(ef5_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango5",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

ef6_95 <- ef6[ef6 >= cuantil_95]

hist(ef6_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango6",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

ef7_95 <- ef7[ef7 >= cuantil_95]

hist(ef7_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango7",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

ef8_95 <- ef8[ef8 >= cuantil_95]

hist(ef8_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango8",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

ef9_95 <- ef9[ef9 >= cuantil_95]

hist(ef9_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango9",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

ef10_95 <- ef10[ef10 >= cuantil_95]

hist(ef10_95,
     main = "Valores del Cuantil 95 del Estadístico F Rango10",
     xlab = "Valores del Estadístico F",
     ylab = "Frecuencia",
     col = "blue",
     border = "black",
     breaks = 50)

curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
      col = "darkgreen", 
      lwd = 2, 
      add = TRUE, 
      from = cuantil_95)

# Grafica solo cuantil 95

par(mfrow = c(2, 5))
for (i in 1:10) {
  ef1_95 <- get(paste0("ef",i,"_95"))  
  hist(ef1_95,
       main = paste("ef",i),
       xlab = "Valores del Estadístico F",
       ylab = "Frecuencia",
       col = "blue",
       border = "black",
       breaks = 50)
  
  curve(df(x, df1 = gl_numerador, df2 = gl_denominador), 
        col = "darkgreen", 
        lwd = 2, 
        add = TRUE, 
        from = cuantil_95)
}

par(mfrow = c(1, 1))

Distancia de Hellinger

d1 = HellingerDist(Fd(df1 = gl_numerador , df2 = gl_denominador),e2 = ef1)
d2 = HellingerDist(Fd(df1 = gl_numerador , df2 = gl_denominador),e2 = ef2)
d3 = HellingerDist(Fd(df1 = gl_numerador , df2 = gl_denominador),e2 = ef3)
d4 = HellingerDist(Fd(df1 = gl_numerador , df2 = gl_denominador),e2 = ef4)
d5 = HellingerDist(Fd(df1 = gl_numerador , df2 = gl_denominador),e2 = ef5)
d6 = HellingerDist(Fd(df1 = gl_numerador , df2 = gl_denominador),e2 = ef6)
d7 = HellingerDist(Fd(df1 = gl_numerador , df2 = gl_denominador),e2 = ef7)
d8 = HellingerDist(Fd(df1 = gl_numerador , df2 = gl_denominador),e2 = ef8)
d9 = HellingerDist(Fd(df1 = gl_numerador , df2 = gl_denominador),e2 = ef9)
d10 = HellingerDist(Fd(df1 = gl_numerador, df2 = gl_denominador),e2 = ef10)
tabladistancias = as.data.frame(rbind(d1,d2,d3,d4,d5,d6,d7,d8,d9,d10))

barplot(tabladistancias$`Hellinger distance`)