Pruebas estadísticas para validez de los geneadores de números Pseudoaleatorios (GNP)

Para la validez estadística de los GNP, vamos a generar 50 números Pseudoaleatorios con el método congruencial lineal mixto. probaremos la aleatoriedad mediante la prueba \(\chi^2\) y la prueba de kolmogorov-Smirnoff , y para la prueba de independencia estadística aplicaremos las pruebas de las Corridas.

Generando los 50 números U_i Pseudoaleatorios con el método secuencial mixto \(X_{i+1}= aX_i+c\,\, mod\,\, m\)

a <- 317 
c <- 15  
m <- 571
X_n <- 41 # semilla
random.number<-numeric(50) # vector numérico de longitud 50
 for (i in 1:50)
   {X_n<-(a*X_n+c)%%m
   random.number[i]<-X_n/m # números en el intervalo [0,1]
   }
 random.number
##  [1] 0.78809107 0.85113835 0.83712785 0.39579685 0.49387040 0.58318739
##  [7] 0.89667250 0.27145359 0.07705779 0.45359019 0.81436077 0.17863398
## [13] 0.65323993 0.10332750 0.78108581 0.63047285 0.88616462 0.94045534
## [19] 0.15061296 0.77057793 0.29947461 0.95971979 0.25744308 0.63572680
## [25] 0.55166375 0.90367776 0.49211909 0.02802102 0.90893170 0.15761821
## [31] 0.99124343 0.25043783 0.41506130 0.60070053 0.44833625 0.14886165
## [37] 0.21541156 0.31173380 0.84588441 0.17162872 0.43257443 0.15236427
## [43] 0.32574431 0.28721541 0.07355517 0.34325744 0.83887916 0.95096322
## [49] 0.48161121 0.69702277
  1. Pruebas de Aleatoriedad:

1.1 Prueba \(\chi^2\)

Se divide el rango de los números aleatorios en \(k\) subintervalos y se cuenta cuántos números caen en cada intervalo. se compara las frecuencias observadas \(O_i\) con las frecuencias esperadas \(E_i\) usando la distribución Chi-cuadrado. Cálculo de la estadística \(\chi^2\) \[ \chi^2 = \sum_{i=1}^{k} \frac{(O_i - E_i)^2}{E_i} \]

Si \(\chi^2 < \chi_{(\alpha,k-1)}\), no se rechaza que el conjunto de números Pseudoaleatorio provienen de una distribución uniforme [0,1].

hist(random.number, breaks=6)$count # obteniendo las frecuencias observadas con 6 intervalos
## [1] 10 10  9  8 13
 library(agricolae) # cargando la libreria agricolae 
## Warning: package 'agricolae' was built under R version 4.3.3
 histo <- hist(random.number, breaks=6)

 Tabla <- table.freq(histo) 
 lim_inf <-Tabla$Lower #limite inferior del intervalo
 lim_sup <- Tabla$Upper # limite superior del intervalo
 
 obser <- Tabla$Frequency
 Ei <- length(random.number)/length(obser) # Valor esperado en una uniforme es E= n/#intervalos
 cbind(lim_inf,lim_sup,obser,Ei) # visualizacion de las frecuencias observadas y esperadas en los intervalos.
##      lim_inf lim_sup obser Ei
## [1,]     0.0     0.2    10 10
## [2,]     0.2     0.4    10 10
## [3,]     0.4     0.6     9 10
## [4,]     0.6     0.8     8 10
## [5,]     0.8     1.0    13 10
 # calculando la estadística chi cuadrada
  x2 <- (obser-Ei)^2/Ei 
 x2
## [1] 0.0 0.0 0.1 0.4 0.9
 chicuad <- sum(x2)
 chicuad # Estadistica chi cuadrada
## [1] 1.4
 dchi <- qchisq(0.05,length(obser)-1,lower.tail=F)# valor de la distribución chi cuadrada con k-1 grado de libertad y nivel de signiicancia 0.05
 
 # Decisión estadística
 
 ifelse(chicuad < dchi,"Los $U_i$ provienen de una distribución uniforme [0,1]", "Los $U_i no siguen una uniforme [0,1] ")
## [1] "Los $U_i$ provienen de una distribución uniforme [0,1]"

1.2 Prueba de Kolmogorov-Smirnov (Prueba K-S).

Se utiliza para comparar si una muestra sigue una distribución específica (normal, exponencial,uniforme etc.).

La hipótesis nula \((H_0 )\) es que la muestra proviene de la distribución teórica especificada, para nuetro caso la distribución teórica , es la uniforme [0,1]

La distribución Teórica: Para una distribución uniforme en el intervalo \([0, 1]\), la función de distribución acumulada es: \[ F(x) = x \] La estadística de prueba: \[ D_n = \max |F_n(x) - F(x)| \]

Si \(D_n < D_{\alpha, n}\) , se concluye que los datos provienen de una distribución uniforme. En R , empleamos la función \(\textit{ks.test}\)

Aplicando la prueba K-S

# Probando si los U_i provienen de una uniforme [0,1]
test_ks <- ks.test(random.number,"punif",0,1)
test_ks
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  random.number
## D = 0.090578, p-value = 0.7728
## alternative hypothesis: two-sided
# Decision
 test_ks$p.value # P válido
## [1] 0.7727792
ifelse(test_ks$p.value < 0.05, " Los u_i no siguen una distribución uniforme","Los u_i siguen la distribución uniforme [0,1]") # 
## [1] "Los u_i siguen la distribución uniforme [0,1]"
  1. Prueba de Independencia. Otro aspecto a validar es la independencia de los números Psuedoaleatorios generados. Pruebas estadísticas que tratan de corroborar si los números en el intervalo [0, 1] son independientes.

Hipótesis:
\(H_0 :\) El conjunto de números Pseudoaleatorios son independientes.
\(H_1 :\) El conjunto de números Pseudoaleatorios NO son independientes

2.1 Prueba de Corrida: arriba y abajo

El procedimiento de esta prueba consiste en determinar una secuencia \(S\) binaria de unos y ceros, de acuerdo con una comparación entre los números aleatorios \(u_i\) y \(u_{i+1}\) Después se determina el número de corridas observadas, \(C_0\) (una corrida se identifica como la cantidad de unos o ceros consecutivos)

El procedimiento de esta prueba consiste en determinar una secuencia de números (S) que sólo contiene unos y ceros, de acuerdo a la siguiente comparación

\[ \left \{ \begin{aligned} 1 & \quad \text{si}\, u_i >u_{i+1}\\ 0 & \quad \text{si}\, u_i \leq u_{i+1} \end{aligned} \right . \] - Determinamos el total de corridas \(C_0\), observadas.

Si \(|Z_0| > Z_{\frac{\alpha}{2}}\), se concluye que los números del conjunto \(u_i\) no son independientes. De lo contrario no se puede rechazar que el conjunto de \(u_i\) sea independiente.

diff(random.number) # halla la diferencia entre un valor  y el otro en un vector
##  [1]  0.06304729 -0.01401051 -0.44133100  0.09807356  0.08931699  0.31348511
##  [7] -0.62521891 -0.19439580  0.37653240  0.36077058 -0.63572680  0.47460595
## [13] -0.54991243  0.67775832 -0.15061296  0.25569177  0.05429072 -0.78984238
## [19]  0.61996497 -0.47110333  0.66024518 -0.70227671  0.37828371 -0.08406305
## [25]  0.35201401 -0.41155867 -0.46409807  0.88091068 -0.75131349  0.83362522
## [31] -0.74080560  0.16462347  0.18563923 -0.15236427 -0.29947461  0.06654991
## [37]  0.09632224  0.53415061 -0.67425569  0.26094571 -0.28021016  0.17338004
## [43] -0.03852890 -0.21366025  0.26970228  0.49562172  0.11208406 -0.46935201
## [49]  0.21541156
S<-ifelse(diff(random.number) > 0, 1, 0) # seuencia de ceros y unos
S
##  [1] 1 0 0 1 1 1 0 0 1 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0 1 0 0 1 0 1 0 1 1 0 0 1 1 1
## [39] 0 1 0 1 0 0 1 1 1 0 1
# devuelve 1 si el número es mayor que el anterior y 0 en caso contrario.

# Detectar cambios
cambios <- abs(diff(S))

# Contar las corridas
corridas <- sum(cambios) + 1
corridas
## [1] 35
mu = (2*length(random.number)-1)/3 # media espearada de corrida 
mu
## [1] 33
varianza <- (16*length(random.number)-29)/90
 
Z <- (corridas-mu)/sqrt(varianza) # valor de la estadística  z
Z
## [1] 0.6833199
# Decisión estadística 

ifelse(Z<1.96,"Los u_i son independientes","Los u_i son dependeientes")
## [1] "Los u_i son independientes"