Search code examples
rplotrandom

Spectral Test in R

I would like to change my code in a way that 10000 points are placed in a [0,1]^2 plot. I tries changing 256 to 10000 but it generates weird placement. I should change the factors 137 and 187 but not sure how to change it. Anyone knows the logic behind?

Working sample:

nSim = 256
X=rep(0,nSim)
for (i in 2:nSim){
    X[i] = (137*X[i-1]+187)%%256 
}
plot(X[-1],X[-nSim],col="blue",type="p",pch="x",lwd=2)

enter image description here

My code:

nSim = 10000
X=rep(0,nSim)
for (i in 2:nSim){
  X[i] = ((137*X[i-1]+187)%%nSim)
}
plot(X[-1]/nSim,X[-nSim]/nSim,col="blue", type="p",pch=20,lwd=2)

enter image description here

Solution

  • What you have there is called a linear congruential generator taking a general form

    Xn+1 = (aXn+c) mod m

    The goal is to choose such a, c, and m that the generated sequence would resemble a random one as much as possible. There is no best way to pick a, c, and m, but there is something we can easily do.

    You already have chosen m = 10000. Then we may use a well-known theorem (see, e.g., the end of this post) how to pick such a and c that the generated numbers would start repeating themselves only after m steps.

    Condition 1: c and m need to be relatively prime. we have that m = 10000 = 2454. Meanwhile, 187 is not divisible by neither 2 nor 4, so we are good.

    Condition 2: if 4 divides m, then 4 also needs to divide a-1. In our case 137-1 = 136 is divisible by 4: 136/4 = 34, so we are good here too.

    Condition 3: if any prime p divides m, then p also needs to divide a-1. We already checked p=2 in the previous step, so we are left with p=5. But 5 does not divide 137-1 = 136! Indeed, as a result of this we get

    length(unique(X))
# [1] 2000

    Namely, in a sequence of length 10000 we have only 2000 unique numbers, meaning that the same numbers are repeated five times. So, then we need such a that a-1 would be divisible by 4 and 5. That gives quite a few options! We may pick, for instance, a = 4*5*6 + 1 = 121.

    Thus, using a = 121, c = 187, and m = 10000 gives

    length(unique(X))
# [1] 10000

    and a plot

    enter image description here

    It's still somewhat regular, but definitely better than the one before. You may keep experimenting with different a and c that satisfy the three conditions.