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

library(statip)

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

## 
## Attaching package: 'statip'

## The following object is masked from 'package:distrEx':
## 
##     plot

## The following object is masked from 'package:distr':
## 
##     plot

library(tidyr)

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

cc = 1000376863
set.seed(cc)
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
}

print(paste("El tamaño de la grilla es", a, "x", b))

## [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.006190945
## 
## $expected
## [1] -0.006993007
## 
## $sd
## [1] 0.008379577
## 
## $p.value
## [1] 0.9237459

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
}

calcular_moran <- function(tabla) {
  sapply(tabla[, 3:102], function(columna) Moran.I(columna, w)$p.value)
}

# Aplicar la función a cada tabla
resultados_moran <- lapply(tablas, calcular_moran)

# Crear un data frame con los p-valores 
df_pvalores <- data.frame(
  pvalores = unlist(resultados_moran),
  tabla = rep(1:length(resultados_moran), each = 100) 
)


ggplot(df_pvalores, aes(x = pvalores)) +
  geom_histogram(binwidth = 0.05, color = "black", fill = "skyblue") +
  geom_vline(xintercept = 0.05, color = "red", linetype = "dashed", linewidth = 1) +  
  
  labs(title = "Distribución de p-valores del Índice de Moran por Rango",
       x = "p-valor",
       y = "Frecuencia") +
  theme_minimal() +
  facet_wrap(~ tabla, ncol = 2) 

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)$p.value
m2 <- Moran.I(r2$media, w)$p.value
m3 <- Moran.I(r3$media, w)$p.value
m4 <- Moran.I(r4$media, w)$p.value
m5 <- Moran.I(r5$media, w)$p.value
m6 <- Moran.I(r6$media, w)$p.value
m7 <- Moran.I(r7$media, w)$p.value
m8 <- Moran.I(r8$media, w)$p.value
m9 <- Moran.I(r9$media, w)$p.value
m10 <- Moran.I(r10$media, w)$p.value
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)
abline(h = 0.05, col = "red")

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

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

# 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 
## 1.0359573 1.0069870 0.7140450 0.4941639 1.6153835 0.7647960 0.6475578 0.3776863 
##      sim9     sim10     sim11     sim12     sim13     sim14     sim15     sim16 
## 0.2405995 1.0715126 1.3625834 0.4661630 0.7456730 0.8544826 0.2268002 0.4269961 
##     sim17     sim18     sim19     sim20     sim21     sim22     sim23     sim24 
## 0.8206079 1.1603940 0.3042917 0.5317910 0.8081928 1.5402828 0.5923557 1.3595678 
##     sim25     sim26     sim27     sim28     sim29     sim30     sim31     sim32 
## 0.8528755 0.9160554 0.9523219 0.6770102 2.2647749 1.4649582 0.4764506 0.3240481 
##     sim33     sim34     sim35     sim36     sim37     sim38     sim39     sim40 
## 0.6012303 0.4050221 1.4084562 0.9985730 1.8003595 0.2987230 0.6969862 0.5593294 
##     sim41     sim42     sim43     sim44     sim45     sim46     sim47     sim48 
## 0.4444674 1.3353080 1.0197188 0.7222278 1.2501834 1.4398292 0.4346018 0.9425081 
##     sim49     sim50     sim51     sim52     sim53     sim54     sim55     sim56 
## 0.7145364 1.2075784 0.7033118 0.8339556 0.9048862 1.5557245 1.3319325 1.2460926 
##     sim57     sim58     sim59     sim60     sim61     sim62     sim63     sim64 
## 0.6999305 1.2439526 0.6934259 2.5542595 1.5595637 1.3507918 0.6562236 0.5253461 
##     sim65     sim66     sim67     sim68     sim69     sim70     sim71     sim72 
## 1.6121199 1.1245503 0.4963043 0.7621490 0.7582430 1.2232827 1.1250761 0.9985771 
##     sim73     sim74     sim75     sim76     sim77     sim78     sim79     sim80 
## 1.7063430 0.6465764 0.8619237 0.9066171 0.4895844 1.1380776 0.7134978 1.2306241 
##     sim81     sim82     sim83     sim84     sim85     sim86     sim87     sim88 
## 0.9409849 1.5693282 0.2219053 1.1247629 0.7188291 0.8965573 0.8399943 1.4481618 
##     sim89     sim90     sim91     sim92     sim93     sim94     sim95     sim96 
## 0.7958536 0.8833244 0.5316503 0.5212424 1.1701936 1.0870618 1.2413254 0.7683567 
##     sim97     sim98     sim99    sim100 
## 2.3462536 0.5841306 1.0222248 0.5443690

# 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 PARA DISTRIBUCIÓN COMPLETA

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`)
abline(h = mean(tabladistancias$`Hellinger distance`, col = "red"))

Distancias de Hellinger en cada DECIL

Visualización de la segmentación en Deciles

# Crear una secuencia de valores para la variable x
x_vals <- seq(0, 5, length.out = 100)

# Calcular la densidad de la distribución F para cada valor de x
dens_vals <- df(x_vals, df1 = gl_numerador, df2 = gl_denominador)

# Crear un data frame con los valores de x y su densidad correspondiente
df_plot <- data.frame(x = x_vals, density = dens_vals)

# Calcular los deciles (10%, 20%, ..., 90%)
deciles <- qf(seq(0.1, 1, by = 0.1), df1 = gl_numerador, df2 = gl_denominador)

# Visualizar la distribución F con la segmentación en deciles
ggplot(df_plot, aes(x = x, y = density)) +
  geom_line(color = "blue", linewidth = 1) + # Línea de la función de densidad
  geom_area(data = df_plot %>% filter(x <= deciles[1]), fill = "red", alpha = 0.4) + 
  geom_area(data = df_plot %>% filter(x > deciles[1] & x <= deciles[2]), fill = "orange", alpha = 0.4) + 
  geom_area(data = df_plot %>% filter(x > deciles[2] & x <= deciles[3]), fill = "yellow", alpha = 0.4) + 
  geom_area(data = df_plot %>% filter(x > deciles[3] & x <= deciles[4]), fill = "green", alpha = 0.4) + 
  geom_area(data = df_plot %>% filter(x > deciles[4] & x <= deciles[5]), fill = "lightblue", alpha = 0.4) + 
  geom_area(data = df_plot %>% filter(x > deciles[5] & x <= deciles[6]), fill = "blue", alpha = 0.4) + 
  geom_area(data = df_plot %>% filter(x > deciles[6] & x <= deciles[7]), fill = "purple", alpha = 0.4) + 
  geom_area(data = df_plot %>% filter(x > deciles[7] & x <= deciles[8]), fill = "pink", alpha = 0.4) + 
  geom_area(data = df_plot %>% filter(x > deciles[8] & x <= deciles[9]), fill = "brown", alpha = 0.4) +
  geom_area(data = df_plot %>% filter(x > deciles[9] & x <= deciles[10]), fill = "black", alpha = 0.4) +
  labs(title = "Distribución F segmentada en deciles", 
       x = "Valores de la distribución F", 
       y = "Densidad") +
  theme_minimal()

# Crear la secuencia de valores para la distribución F
x_vals <- seq(0, 5, length.out = 100)
dens_vals <- df(x_vals, df1 = gl_numerador, df2 = gl_denominador)

# Asumimos que las distribuciones ef1, ef2, ..., ef10 ya están en el entorno
ef_list <- list(ef1, ef2, ef3, ef4, ef5, ef6, ef7, ef8, ef9, ef10)

# Calcular los deciles para la distribución F
deciles <- qf(seq(0.1, 0.9, by = 0.1), df1 = gl_numerador, df2 = gl_denominador)

# Función para calcular la distancia de Hellinger en un decil, con manejo de errores
hellinger_by_decil <- function(lower_bound, upper_bound, f_dist, ef_dist) {
  # Filtrar los valores en el rango del decil
  f_decile <- f_dist[x_vals >= lower_bound & x_vals <= upper_bound]
  ef_decile <- ef_dist[x_vals >= lower_bound & x_vals <= upper_bound]
  
  # Normalizar las distribuciones en el rango para compararlas
  f_decile <- f_decile / sum(f_decile)
  ef_decile <- ef_decile / sum(ef_decile)
  
  # Calcular la distancia de Hellinger con manejo de errores
  hellinger_dist <- tryCatch({
    hellinger(f_decile, ef_decile)
  }, error = function(e) {
    warning("Error en el cálculo de Hellinger: ", e$message)
    return(NA)
  })
  
  return(hellinger_dist)
}

# Crear un dataframe para almacenar los resultados de Hellinger
results <- data.frame(ef_number = integer(), decil_number = integer(), hellinger_value = numeric())

# Calcular la distancia de Hellinger para cada decil y cada distribución ef
for (ef_idx in 1:10) {
  ef_dist <- ef_list[[ef_idx]]
  
  for (i in 1:length(deciles)) {
    if (i == 1) {
      lower_bound <- 0
    } else {
      lower_bound <- deciles[i-1]
    }
    upper_bound <- deciles[i]
    
    # Calcular la distancia de Hellinger en este decil
    h_dist <- hellinger_by_decil(lower_bound, upper_bound, dens_vals, ef_dist)
    
    # Almacenar el resultado en el dataframe
    if (!is.na(h_dist)) {
      results <- rbind(results, data.frame(ef_number = ef_idx, decil_number = i, hellinger_value = h_dist))
    }
  }
}

# Graficar los resultados: una gráfica para cada decil con barras y línea de la media
for (i in 1:9) {
  # Filtrar los resultados para el decil i
  decil_data <- results %>% filter(decil_number == i)
  
  # Calcular la media de las distancias de Hellinger en este decil
  mean_hellinger <- mean(decil_data$hellinger_value, na.rm = TRUE)
  
  # Crear la gráfica de barras con la línea de la media
  p <- ggplot(decil_data, aes(x = factor(ef_number), y = hellinger_value)) +
    geom_bar(stat = "identity", fill = "skyblue") +
    geom_hline(yintercept = mean_hellinger, linetype = "dashed", color = "red", linewidth = 1) + # Línea de la media
    scale_x_discrete(labels = 1:10) +  # Escala del eje X de 1 a 10
    labs(title = paste("Distancia de Hellinger - Decil", i),
         x = "Distribución ef", 
         y = "Distancia de Hellinger") +
    theme_minimal()
  
  # Mostrar la gráfica
  print(p)
}