Προσομοίωση

1. Γεννήτριες ψευδοτυχαίων ομοιόμορφων αριθμών - Linear Congruential Generator

Η πιο κάτω συνάρτηση υλοποιεί τον αλγόριθμο Linear Congruential Generator.

lcg <- function(m,a,c,X0, N)
{ # Integer m>1 is the modulus
  # Integer 0<a<m is the multuplier
  # Integer 0<=c<m is the increment
  # X0 natural number or zero is the seed
  # N is the number of uniform samples required
  # returns sequence X1,...,XN of numbers which for good values of m,a,c 
  # behave like uniformly distributed independent random numbers
  x <- c()
  x[1] <- X0
  for(i in 1:N)
  {
    x[i+1] <- (a*x[i] + c)%%m
  }
  return(x[-1])
}

Για παράδειγμα για \(m=8, a=5, c=1, X_0=0, N=10\), παίρνουμε την εξής ακολουθία που συμφωνεί με τον υπολογισμό στη θεωρία.

x <- lcg(8,5,1,0,1000)
x

##    [1] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##   [35] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##   [69] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [103] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [137] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [171] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [205] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [239] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [273] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [307] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [341] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [375] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [409] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [443] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [477] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [511] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [545] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [579] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [613] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [647] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [681] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [715] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [749] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [783] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [817] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [851] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [885] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [919] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [953] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [987] 7 4 5 2 3 0 1 6 7 4 5 2 3 0

y <- lcg(2^32, 1103515245, 12345, 0, 100)
y

##   [1]      12345 3554416254 2802067456 1358247936 2138638336 1459132416
##   [7] 1445521408 2518349824 4044081152 3908327424 3000608768 1742518272
##  [13] 3372910592 4276885504  165046272 4093278272 1942545408 3657516032
##  [19] 3443663872 1679918080 4034749440 4057711616 2286593024 2646606848
##  [25] 3287933952 3973989376  727393280  314999808 2461951040 1193041920
##  [31]  743772160  722993152 2839031808 1267507200 3282063360   60964864
##  [37] 2073391160  706873344 2129500160 4286179328 1536686080  491046912
##  [43]  456742976 2446214016 2724457984 2869986816 2195537408  353867264
##  [49]  784272960 1664139520 1552964864 1405995264 3621569792 2332430336
##  [55] 1786052608 4191236096  984719360 3697364992 1540308992 3166597120
##  [61]   53862400  137400376  200189984 4103570912 3129896960 1959124992
##  [67] 1115721728 3701309440 3960451072 1783504896 3506688000  134201344
##  [73] 3563384896 1101921280 3789323264 1084816384 3971709952 3550215168
##  [79] 2285009920 1195834368 3791748096 2523528192 1038646272  432339968
##  [85] 3852958784  699111424 2522867712 2944792576 1570424832 3017234432
##  [91] 1065396224  754432000 4018767872  352262144  971860032 3780711168
##  [97] 3367454720 4261009408 1395730432 1444820992

Για να πάρουμε δείγμα από την ομοιόμορφη συνεχή στο διάστημα \([0,1]\), επιστρέφουμε \(U_n=\frac{X_n+1}{m+1}\).

lcg01 <- function(m, a, c, X0, N)
{ # Integer m>1 is the modulus
  # Integer 0<a<m is the multiplier
  # Integer 0<=c<m is the increment
  # X0 natural number or zero is the seed
  # N is the number of uniform samples required
  # returns sequence X1,...,XN of numbers which for good values of m,a,c 
  # behave like uniformly distributed independent random numbers
  x <- c()
  x[1] <- X0
  for(i in 1:N)
  {
    x[i+1] <- (a*x[i] + c) %% m
  }
  return((x[-1] + 1) / (m + 1))
}

x <- lcg01(8,5,1,0,1000)
x[1:30]

##  [1] 0.2222222 0.7777778 0.8888889 0.5555556 0.6666667 0.3333333 0.4444444
##  [8] 0.1111111 0.2222222 0.7777778 0.8888889 0.5555556 0.6666667 0.3333333
## [15] 0.4444444 0.1111111 0.2222222 0.7777778 0.8888889 0.5555556 0.6666667
## [22] 0.3333333 0.4444444 0.1111111 0.2222222 0.7777778 0.8888889 0.5555556
## [29] 0.6666667 0.3333333

hist(x)

plot.ts(x, main="LCG Path m=8, a=5, c=1, x0=0")