Η πιο κάτω συνάρτηση υλοποιεί τον αλγόριθμο 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")