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 congruencial 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]
   }
 round(random.number,5)
##  [1] 0.78809 0.85114 0.83713 0.39580 0.49387 0.58319 0.89667 0.27145 0.07706
## [10] 0.45359 0.81436 0.17863 0.65324 0.10333 0.78109 0.63047 0.88616 0.94046
## [19] 0.15061 0.77058 0.29947 0.95972 0.25744 0.63573 0.55166 0.90368 0.49212
## [28] 0.02802 0.90893 0.15762 0.99124 0.25044 0.41506 0.60070 0.44834 0.14886
## [37] 0.21541 0.31173 0.84588 0.17163 0.43257 0.15236 0.32574 0.28722 0.07356
## [46] 0.34326 0.83888 0.95096 0.48161 0.69702
  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"

Haciedo Uso de la librería DescTolls y la función Runstest

library(DescTools)
## Warning: package 'DescTools' was built under R version 4.3.3
RunsTest(S)
## 
##  Runs Test for Randomness
## 
## data:  S
## z = 2.701, runs = 35, m = 22, n = 27, p-value = 0.006913
## alternative hypothesis: true number of runs is not equal the expected number

2.2 Prueba Poker

Los números pseudoaleatorios \(u_i\) se agrupan en combinaciones de varios dígitos decimales (típicamente grupos de 5 decimales). Los grupos de dígitos se clasifican de manera similar a las manos en el juego de póker:

  1. Todo diferente (TD): Ningún decimal se repite en el grupo (Ej: 0.12345).
    la probabilidad de este evento es: \[P(TD)=0.3024\]
  2. Un par (1p): Un decimal aparece dos veces y los demás son diferentes (Ej: 0.14345).

\[P(1p)=0.5040\]

  1. Dos pares (2p): Dos decimales aparecen dos veces cada uno (Ej: 0.14224).
    \[P(2p)=0.1080\]

  2. Tercia (T): Un decimal aparece tres veces, y los otros dos son diferentes (Ej: 0.41443).
    \[P(T)=0.0720\]

  3. Full House (FH): Un dígito aparece tres veces y otro dígito dos veces (Ej: 0.32223).
    \[P(FH)=0.0090\]

  4. Poker (p): Un dígito aparece cuatro veces (Ej: 0.55515).
    \[P(p)=0.0045\]

  5. Quintilla: Todos los dígitos son iguales (Ej: 0.11111).
    \[P(Q)=0.0001\]

library(randtoolbox)
## Loading required package: rngWELL
## This is randtoolbox. For an overview, type 'help("randtoolbox")'.
poker.test(random.number,nbcard=5)
## 
##           Poker test
## 
## chisq stat = 3.9, df = 4, p-value = 0.43
## 
##       (sample size : 50)
## 
##  observed number  0 0 7 2 1 
##  expected number  0.016 0.96 4.8 3.8 0.38