Metodo MIDDLE - SQUARE di Von Neumann

L’idea è di considerare un numero iniziale, detto SEED, indicato con X; X è composto di N cifre e viene scelto in modo arbitrario.

Eleva X al quadrato e considera come nuova simulazione le sue N cifre centrali; poi ripeti a partire dal numero così simulato.

Il limite del metodo sta nel fatto che è molto sensibile alla scelta iniziale di (X,N) e che non è in grado di simulare nuovi numeri se viene simulato uno 0 (00, 000, …).

# N = 2
# X = 81
X = 81
random.numbers = c()
repeat {
  output <- X^2
  u <- substring(as.character(output),
                 ceiling(nchar(as.character(X^2))/2),
                 ceiling(nchar(as.character(X^2))/2)+1)
  X <- as.numeric(u)
  random.numbers = c(random.numbers, X)
  if (length(random.numbers)>9) break
}
random.numbers
##  [1] 56 13 69 76 77 92 46 11 21 41

Metodo Congruenziale LINEARE

Fissato un punto di partenza x ed un divisore m, il generatore simula numeri pseudo-casuali y = x mod m (i.e. y è il resto della divisione x/m).

# genero 2 numeri pseudo casuali; serve cambiare di volta in volta x0
x0 <- 13 
m <- 15
x0%%m
## [1] 13
x0.new <- 13 + 12
x0.new%%m
## [1] 10

Metodo Congruenziale MISTO (Lehmar)

x[i +1] = (a*x[i] + c) mod m

Si aggiorna il numero generato premoltiplicandolo per a ed aggiungendo la costante c.  Poi si esegue di nuovo l’operazione di modulo e si ottiene in output il resto della divisione.

a <- 41
c <- 17
x0 <- 1
m <- 16
x <- c()
Lehmar <- function(x, a, c, m) {
  (a*x + c)%%m
}
x[1] <- Lehmar(x=x0,a,c,m)
for (i in 2:20) {
  x[i] <- Lehmar(x=x[i-1],a,c,m)
}
# nota che il periodo della sequenza è pari ad m (dal 16esimo elemento in poi, la sequenza si ripete)

# per ottenere numeri uniformi, semplicemente divido la sequenza generata per il massimo numero generabile (m) 
u <- x/m; u
##  [1] 0.6250 0.6875 0.2500 0.3125 0.8750 0.9375 0.5000 0.5625 0.1250 0.1875
## [11] 0.7500 0.8125 0.3750 0.4375 0.0000 0.0625 0.6250 0.6875 0.2500 0.3125
#qual è il ciclo della sequenza generata con: a=41, m=32768, c=901, x0=1?
x[1] <- Lehmar(x=1, a=41, m=32768, c=901)
for (i in 2:32770) {
  x[i] <- Lehmar(x=x[i-1], a=41, m=32768, c=901)
}
which(duplicated(x)==1)
## [1] 32769 32770
# risposta: m!

Metodo Congruenziale MOLTIPLICATIVO

Meno efficiente rispetto al metodo Lehmar, consiste nel generare numeri con la formula ricorsiva: x[i+1] = (a*x[i]) mod m

Poi si divide la sequenza generata per m, in modo da ottenere numeri pseudo-casuali uniformi.

m = 32768
a = 61
x0 = 1
x = c()
x[1] <- (a*x0)%%m
for (i in 2:20) {
  x[i] <- (a*x[i-1])%%m
}
u <- x/m; u
##  [1] 0.001861572 0.113555908 0.926910400 0.541534424 0.033599854 0.049591064
##  [7] 0.025054932 0.528350830 0.229400635 0.993438721 0.599761963 0.585479736
## [13] 0.714263916 0.570098877 0.776031494 0.337921143 0.613189697 0.404571533
## [19] 0.678863525 0.410675049

Metodo Congruenziale ADDITIVO

Stavolta il numero x generato non dipende solo dal suo predecessore, ma anche da uno o più dei suoi precedenti ad un passo k, secondo la relazione: x[i+1] = x[i] + x[i-k] mod m

m = 32678
k = 2
x0 = 1
x1 = x0%%m
x2 = x1 + x0%%m 
x <- c()
# serie di fibonacci
x[1] <- x1; x[2] <- x2
for (i in 3:20) {
  x[i] <- x[i-1] + x[i-k]%%m
}
x
##  [1]     1     2     3     5     8    13    21    34    55    89   144   233
## [13]   377   610   987  1597  2584  4181  6765 10946

Test di Casualità

Dopo aver prodotto la sequenza di numeri pseudo-casuali nell’intervallo continuo (0,1), dobbiamo verificare che si tratta davvero di numeri casuali, i.e. numeri INDIPENDENTI estratti da una distribuzione U(0,1).

Il test di casualità si compone quindi di due parti:

  1. Uniformità
  2. Indipendenza

Test di Uniformità

Possiamo performarlo in almeno 3 modi: empirico, Chi-quadrato, KS

A. Test Empirico

(i)  confronto momenti campionari e teorici
(ii) confronto grafico 
set.seed(123)
n <- 10000
u <- runif(n)
u1 <- u[-n]
u2 <- u[-1]
# momenti teorici di U(0,1)
c(1/2, 1/12)
## [1] 0.50000000 0.08333333
# momenti empirici di u 
c(mean(u), var(u))
## [1] 0.49754937 0.08219329
# u si dispone come un quadrato su R2 ?
plot(x=u1, y=u2)

# la densità di u sembra uniforme?
hist(plot=TRUE, x=u, freq=F)
abline(h=1, col='red')

B. Test del Chi-Quadrato

H0: i numeri (u[1],…,u[n]) sono estratti da una pdf Uniforme(0,1)

Confronto le frequenze teoriche con quelle osservate con la stat. test Chi-Quadrato (media degli scostamenti al quadrato).

# divido il supporto (0,1) in k=10 sottointervalli
k = 10
istogramma = hist(x=u, breaks=k, plot=F)
# calcolo frequenze osservate e attese 
freq.obs = istogramma$counts
freq.exp = rep(n/k, k)
# calcolo statistica test Chi-Quadro: k/N * sum((scostamenti)^2)
Chi.squared = k/n*sum((freq.obs - freq.exp)^2)
# confronto con il quantile a livello 95% (o calcolo p-value)
Chi.squared > qchisq(.95, k-1)
## [1] FALSE
pvalue <- pchisq(Chi.squared, df=k-1, lower.tail = F)
pvalue
## [1] 0.2222474

C. Test di KOLMOGOROV - SMIRNOV

Calcolo la max.distanza assoluta d.n tra la fdr teorica e quella empirica. La fdr teorica è, sotto H0, un’uniforme (0,1). La fdr empirica è pari a i/n con i=0,1,..,n per ogni x(i) (valori x riordinati in modo crescente). Sfrutto un risultato asintotico (n->Inf) per calcolare il valore critico della stat. test sqrt(n)*d.n

# riordino i valori di u
x.H0 <- sort(u, decreasing=F)

# definisco fdr empirica 
efdr <- knots(ecdf(x=u)); plot(x.H0,efdr)

# definisco fdr teorica 
tfdr <- 1:n/n

# calcolo statistica test 
d.n <- max(abs(efdr - tfdr))
stat.test <- sqrt(n)*d.n

# calcolo p-value dal risultato asintotico
p.value <- 2*sum((-1)^(1:n-1)*exp(-2*(1:n)^2*stat.test^2))
p.value
## [1] 0.6453502
# alto p.value -> non rifiuto H0: dati uniformi

# funzione R
ks.test(x=u, y=punif)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  u
## D = 0.0073921, p-value = 0.6454
## alternative hypothesis: two-sided

Possiamo concludere dicendo che il comando R runif(n) genera davvero un campione di n osservazioni uniformi e indipendenti: ipotesi confermata da tutti e 3 tipi di test di casualità.