Simulación Muestras sistemáticas

Muestreo Estadístico

Author

Affiliation

Andrés Fabián Rondón Vargas

Departamento de Estadística, Universidad Nacional de Colombia

---

Función de Varianza

\[ Var(\hat{t_\pi}) = a\sum_{\xi = 0}^{a} \left({t}_{m\xi} - \bar{{t}} \right)^2 \]

varianza_pi <- function(muestras_df) {
  a <- ncol(muestras_df)
  # Suma total de cada columna
  totales <- sapply(muestras_df, function(col) sum(col, na.rm = TRUE))
  # Promedio de los totales
  promedio_total <- mean(totales)
  # Varianza pi
  var_pi <- a * sum((totales - promedio_total)^2)
  return(var_pi)
}

Función SCE

\[ \text{SCE} = \sum_{j=1}^{a} n_{\xi} \left( \bar{y}_{\xi} - \bar{y_U} \right)^2 \]

SCE <- function(muestras_df) {
  YbarraU <- mean(as.matrix(muestras_df), na.rm = TRUE) # Media Global
  # Inicializar SCE
  SCE_total <- 0
  # Recorrer cada columna (muestra sistemática)
  for (j in seq_along(muestras_df)) {
    muestra <- muestras_df[[j]]
    n_j <- sum(!is.na(muestra))  # número de elementos válidos
    Ybarra_j <- mean(muestra, na.rm = TRUE) # promedio de la columna
    SCE_total <- SCE_total + n_j * (Ybarra_j - YbarraU)^2
  }
  
  return(SCE_total)
}

Función muestra sistemática.

muestras_sist <- function(poblacion, a) {
  N <- length(poblacion)
  
  # Calcular tamaño base de cada columna
  base_tam <- floor(N / a)
  extras <- N %% a  # Cuántas columnas deben tener una fila más
  
  muestras <- vector("list", a)
  idx <- 1  # índice actual en la población
  
  for (j in 1:a) {
    tam_j <- base_tam + ifelse(j <= extras, 1, 0)  # tamaño de esta columna
    muestras[[j]] <- poblacion[idx:(idx + tam_j - 1)]
    idx <- idx + tam_j
  }
  
  # Convertir a data.frame rellenando con NA para igualar filas
  max_len <- max(sapply(muestras, length))
  muestras_df <- as.data.frame(lapply(muestras, function(x) {
    length(x) <- max_len
    return(x)
  }))
  
  colnames(muestras_df) <- paste0("Muestra_", 1:a)
  return(muestras_df)
}

Simulación

####a) Cuando n es múltiplo de N ####

set.seed(123)
N <- 1000      # Tamaño población
n <- 200       # Tamaño muestra
poblacion <- rnorm(N, mean = 10, sd = 3)

CASO 1 cuando N es múltiplo de n

c <- N %% n
print(paste("Residuo c = ", c))

[1] "Residuo c =  0"

a<-floor(N/n);a

[1] 5

muestras <- muestras_sist(poblacion, a)

Tabla 1: de la muestra sistemática

cuando N es multiplo de n

Muestra_1	Muestra_2	Muestra_3	Muestra_4	Muestra_5
8.318573	16.596431	9.779332	13.222037	11.068850
9.309468	13.937239	6.494046	9.917959	8.025969
14.676125	9.204565	8.095755	9.900009	12.565607
10.211525	11.629582	9.913475	5.451797	13.458809
10.387863	8.756980	12.012088	12.371156	10.828824
15.145195	8.571259	5.048360	9.367797	10.432314

g<-SCE(muestras);g

[1] 6.155054

varianza_pi(muestras)

[1] 6155.054

Ahora vamos a verificar que cuando N es múltiplo de n se cumple que:

\[ Var(\hat{t_\pi})= SCE*N \]

g*N

[1] 6155.054

CASO 2 cuando N es no es múltiplo de n

N <- 21      # Tamaño población
n <- 4       # Tamaño muestra
poblacion <- rnorm(N, mean = 10, sd = 3)

c <- N %% n
print(paste("Residuo c = ", c))

[1] "Residuo c =  1"

a<-floor(N/n);a

[1] 5

muestras <- muestras_sist(poblacion, a)

Tabla 2: de la muestra sistemática

cuando N es distinto de n

Muestra_1	Muestra_2	Muestra_3	Muestra_4	Muestra_5
7.012604	13.12172	8.659122	11.407096	6.214299
6.880135	10.74918	18.392173	9.366259	10.856769
9.946059	17.24862	18.496678	10.561153	15.247742
9.603475	12.05559	6.343865	10.682628	9.507730
2.351972	NA	NA	NA	NA

g<-SCE(muestras);g

[1] 110.665

varianza_pi(muestras)

[1] 1090.134

Ahora vamos a verificar que cuando N no es múltiplo de n no se cumple que:

\[ Var(\hat{t_\pi})= SCE*N \]

g*N

[1] 2323.965

Simulación mas general

n_que_cumplen <- function(N, tolerancia = 1e-8) {
  # Generar población simulada
  set.seed(123)
  poblacion <- rnorm(N, mean = 10, sd = 3)
  
  # Encontrar divisores de N
  divisores <- which(N %% 1:(N - 1) == 0)
  
  # Inicializar vector de resultados
  cumplen <- c()
  
  for (n in 2:(N - 1)) {
    muestras <- muestras_sist(poblacion, n)
    sce <- SCE(muestras)
    varpi <- varianza_pi(muestras)
    
    if (abs(sce * N - varpi) < tolerancia) {
      cumplen <- c(cumplen, n)
    }
  }
  
  # Imprimir resultados
  cat("Los divisores de", N, "son:", divisores, "\n")
  cat("Los valores de n que cumplen que SCE * N = varianza_pi son:", cumplen, "\n")
  
  return(cumplen)
}

set.seed(123)
n_que_cumplen(32)

Los divisores de 32 son: 1 2 4 8 16 
Los valores de n que cumplen que SCE * N = varianza_pi son: 2 4 8 16 

[1]  2  4  8 16

CASO 3: variable auxiliar con N múltiplo de n

# Función que simula varianzas para un N dado (orden ASCENDENTE)
simular_varianzas <- function(N, n = 10) {
  poblacion <- rnorm(N, mean = 10, sd = 3)
  x_aux <- poblacion + rnorm(N, mean = 0, sd = 0.2)  # Alta correlación (~0.98)

  # Sin ordenar
  muestras_no_ord <- muestras_sist(poblacion, n)
  var_no <- varianza_pi(muestras_no_ord)

  # Ordenando por la auxiliar (ASCENDENTE)
  orden <- order(x_aux)  # Cambio clave: orden ascendente
  muestras_ord <- muestras_sist(poblacion[orden], n)
  var_ord <- varianza_pi(muestras_ord)

  return(c(no = var_no, ord = var_ord))
}

# Simular para distintos N
Ns <- seq(50, 1000, by = 50)
resultados <- sapply(Ns, simular_varianzas)

# Graficar
plot(Ns, resultados["no", ], type = "l", col = "red", lwd = 2,
     ylim = range(resultados), ylab = "Varianza del estimador",
     xlab = "Tamaño de la población (N)", main = "Varianza: orden ascendente vs sin ordenar")
lines(Ns, resultados["ord", ], col = "blue", lwd = 2)
legend("topright", legend = c("Ordenado", "Sin Ordenar"),
       col = c("red", "blue"), lwd = 2)