Generador de números aleatorios con congruencias lineales

X(i+1) = a * X(i) + c (mod m)

a: multiplicador

c: incremento

m: módulo

n: número de aleatorios

gnacl <- function(a, c, m, x0, n) {
    x <- numeric(n)
    x[1] = x0
    for (i in 1:n) {
        x[i + 1] = (a * x[i] + c)%%m
    }
    y <- x/m
    return(y)
}
y <- gnacl(20, 0, 53, 21, 60)
y

##  [1] 0.39623 0.92453 0.49057 0.81132 0.22642 0.52830 0.56604 0.32075
##  [9] 0.41509 0.30189 0.03774 0.75472 0.09434 0.88679 0.73585 0.71698
## [17] 0.33962 0.79245 0.84906 0.98113 0.62264 0.45283 0.05660 0.13208
## [25] 0.64151 0.83019 0.60377 0.07547 0.50943 0.18868 0.77358 0.47170
## [33] 0.43396 0.67925 0.58491 0.69811 0.96226 0.24528 0.90566 0.11321
## [41] 0.26415 0.28302 0.66038 0.20755 0.15094 0.01887 0.37736 0.54717
## [49] 0.94340 0.86792 0.35849 0.16981 0.39623 0.92453 0.49057 0.81132
## [57] 0.22642 0.52830 0.56604 0.32075 0.41509

unique(y)

##  [1] 0.39623 0.92453 0.49057 0.81132 0.22642 0.52830 0.56604 0.32075
##  [9] 0.41509 0.30189 0.03774 0.75472 0.09434 0.88679 0.73585 0.71698
## [17] 0.33962 0.79245 0.84906 0.98113 0.62264 0.45283 0.05660 0.13208
## [25] 0.64151 0.83019 0.60377 0.07547 0.50943 0.18868 0.77358 0.47170
## [33] 0.43396 0.67925 0.58491 0.69811 0.96226 0.24528 0.90566 0.11321
## [41] 0.26415 0.28302 0.66038 0.20755 0.15094 0.01887 0.37736 0.54717
## [49] 0.94340 0.86792 0.35849 0.16981

length(unique(y))

## [1] 52

b <- 1:10/10
hist(y, col = "purple")

y <- gnacl(1093, 0, 86436, 12, 2000)
length(unique(y))

## [1] 343

k <- 0:10/10
hist(y, breaks = k, col = "aquamarine1")

plot(y[1:(2000 - 1)], y[2:2000])

length(unique(y))

## [1] 343

y <- gnacl(33, 12, 64, 57, 100)
k <- 0:10/10
hist(y, breaks = k, col = "greenyellow")

plot(y[1:(100 - 1)], y[2:100])

length(unique(y))

## [1] 16

date()

## [1] "Mon Oct  7 19:59:01 2013"