9  Random Number Generation

A sequence of random numbers \(R_1, R_2, \ldots\) must be independent and uniform.

That is each random number \(R\) must be an independent sample drawn from a continous uniform distribution between 0 and 1. Thus, the pdf is

\[ f(x) = \begin{cases} 1 & 0 \leq x \leq 1 \\ 0 & \text{otherwise} \end{cases} \]

The expected value of \(R\) is

\[ E[R] = \int_0^1 x dx = \frac{1}{2} \]

The variance of \(R\) is

\[ Var[R] = E[R^2] - (E[R])^2 = \frac{1}{3} - \frac{1}{4} = \frac{1}{12} \]

9.1 Pseudo-Random Number Generator (PRNG)

A PRNG generates a random sequence starting from a seed value in a certain range and has a period after which the sequence repeats.

A good prng has a large period size, high density, fast generation, and closely approxiamtaes the uniform distribution.

9.2 Linear Congruential Generator (LCG)

This method produces a sequence of random integers between 0 and \(m-1\) using the following recursion:

\[ X_{i+1} = (aX_i + c) \mod m \]

and to convert this to a random number between 0 and 1, we can use

\[ R_i = \frac{X_i}{m} \]

\(X_0\) is the seed value, \(a\) is the multiplier, \(c\) is the increment, and \(m\) is the modulus.

This version is called mixed LCG. Another version is the multiplicative LCG where \(c=0\). That is

\[ X_{i+1} = (aX_i) \mod m \]

9.2.1 Period of the LCG

Because the value of \(X_i\) is a function of \(X_{i-1}\), whenever \(X_i\) repeats, the sequence will start repeating. The length of the cycle until repeating a value is called the period of the LCG and the maximum possible period is \(m\). If the period is \(m\), then the LCG is called a full period LCG.

9.2.2 Some Rules Regarding Maximum Period

9.2.2.1 Case 1

If

  • \(m = 2^b\)
  • \(c \neq 0\) (that is, mixed LCG)
  • \(c\) is relatively prime to \(m\) (that is, \(gcd(c, m) = 1\); in this case it means that \(c\) is odd)
  • \(a = 1 + 4k\) for some integer \(k\) (that is, \(a \mod 4 = 1\))

Then the longest possible period is \(m\).

9.2.2.2 Case 2

If

  • \(m = 2^b\)
  • \(c = 0\) (that is, multiplicative LCG)
  • \(X_0\) is odd
  • \(a = 3 + 8k\) or \(a = 5 + 8k\) for some integer \(k\) (that is, \(a \mod 8 = 3\) or \(a \mod 8 = 5\))

Then the longest possible period is \(m/4\).

9.2.2.3 Case 3

If

  • \(m\) is prime

Then the longest possible period is \(m-1\).

Note

This is what the slides state, but the textbook explains it in more detail.

For \(m\) a prime number and \(c=0\), the longest possible period is \(P=m-1\), which is achieved whenever the multipllier \(a\) has the property that the smallest integer \(k\) such that \(a^k - 1\) is divisible by \(m\) is \(k=m-1\).

The condition on \(a\) can also be stated in a different way: \(a\) is a primitive root modulo \(m\). To check if \(a\) is a primitive root modulo \(m\), we can check if \(a^{(m-1)/p} \mod m \neq 1\) for every prime factor \(p\) of \(m-1\). Here \(m-1\) is the Euler totient function \(\phi(m)\), which counts the number of integers between 1 and \(m\) that are coprime to \(m\).

Another side-info if that in this case, the period is \(m-1\) because the sequence will go through all the integers from 1 to \(m-1\) before repeating, with 0 never appearing in the sequence.

Needless to say, like any multiplicative LCG, \(X_0 \neq 0\) is required to avoid the sequence being stuck at 0.

If the problem does not satisfy any of the above cases, you have to calculate it by looping yourself.

9.3 Combined Linear Congruential Generator

Warning

This section is copied from the textbook.

In the books reasoning for using \((-1)^{j-1}\), there is a mistake in the given example where he subtracted two values rather than add them.

This is atechnique to combine multiple LCGs to get one with a longer period.

They are based on the following theorem:

If \(W_{i,1}, W_{i,2}, \ldots, W_{i,k}\) are any independent, discrete random variables (not necessarily identically distributed), but one of them, say \(W_{i,1}\), is uniformly distributed on the integers from 0 to \(m_j-2\), then

\[ W_i = \left( \sum_{j=1}^k W_{i,j} \right) \mod (m_1 - 1) \]

is uniformly distributed on the integers from 0 to \(m_1 - 2\).

We can apply this results to form a combined generator.

Let \(X_{i,1}, X_{i,2}, \ldots, X_{i,k}\) be the \(i^{th}\) outputs from \(k\) different multiplicative congruential generators, where the \(j^{th}\) generator has prome modulud \(m_j\) and the multipler \(a_j\) is chosen so that the period is \(m_j - 1\). Then, the \(j^{th}\) generator is producing integers \(X_{i,j}\) that are approximately uniformly distributed on the integers from 1 to \(m_j - 1\), and \(W_{i,j} = X_{i,j} - 1\) is uniformly distributed on the integers from 0 to \(m_j - 2\). L’Ecuyer therefore suggests combined generators of the form

\[ X_i = \left( \sum_{j=1}^k (-1)^{j-1} X_{i,j} \right) \mod (m_1 - 1) \]

with

\[ R_i = \begin{cases} \frac{X_i}{m_1} & \text{if } X_i > 0 \\ \frac{m_1 - 1}{m_1} & \text{if } X_i = 0 \end{cases} \]

The maximum possible period of this combined generator is

\[ P = \frac{(m_1 - 1)(m_2 - 1) \cdots (m_k - 1)}{2^{k-1}} \]

Tip

I think \(P\) is the LCM of the periods of the individual generators.

Notice that the period can be larger than \(m_1 - 1\) which we use as the modulus.

Since for the combined generator \(X_i\) isn’t a function of \(X_{i-1}\), repeating a single value doesn’t imply repeating the whole series. Thus, the period can indeed be larger than \(m_1-1\).

Important

The slides make changes from the book!!!

It states that the \(j^{th}\) generator produces that are uniformly distributed on the integers from 0 to \(m_j - 2\), rather than from 1 to \(m_j - 1\) as stated in the book.

Moreover, it sets \(R_i = \frac{m_1+1}{m_1}\) when \(X_i = 0\) which is just ridiculous (a higher than 1 probability).

9.4 Random Streams

You can create multiple random streams using python

import numpy as np

seed_seq = np.random.SeedSequence(12345)
streams = [np.random.default_rng(s) for s in seed_seq.spawn(5)]

for stream in streams:
    print(stream.random(5))

9.5 Uniformity Tests

We want to verify two properties of the generated random numbers:

  1. Uniformity
  2. Independence

The uniformity tests ensure that the generated random numbers approximate the uniform distribution.

We will focus on two tests:

  1. Kolmogorov-Smirnov Test
  2. Chi-Square Test

Notice that both are acceptable for testing the uniformity of a sample of data provided the sample size is large. However, the Kolmogorov-Smirnov test is more powerful and is recommended. Furthermore, it can be applied to small sample sizes, while the chi-square test is valid only for large sample sizes (usually \(N > 50\)).

9.5.1 Kolmogorov-Smirnov Test

  1. Generate \(n\) random numbers and sort them in increasing order: \(R_{(1)} \leq R_{(2)} \leq \ldots \leq R_{(n)}\).
  2. Compute \[ D^+ = \max_{1 \leq i \leq n} \left( \frac{i}{n} - R_{(i)} \right) \] \[ D^- = \max_{1 \leq i \leq n} \left( R_{(i)} - \frac{i-1}{n} \right) \]
  3. Compute \(D = \max(D^+, D^-)\).
  4. Locate the critical value \(D_\alpha\) from the Kolmogorov distribution table for the desired significance level \(\alpha\) and sample size \(n\).
  5. If \(D > D_\alpha\), reject the null hypothesis that the numbers are uniformly distributed. Otherwise, fail to reject the null hypothesis.

9.5.2 Chi-Square Test

  1. Divide the interval [0, 1) into \(n\) equal classes. Thus, \(E_i = \frac{N}{n}\) for each class, where \(N\) is the total number of generated random numbers.
  2. For \(N\) observations, count the number of already obserced observations in each class: \(O_i\).
  3. Compute the chi-square statistic: \[ \chi^2 = \sum_{i=1}^n \frac{(O_i - E_i)^2}{E_i} \]
  4. Locate the critical value \(\chi^2_{\alpha,n-1}\) from the chi-square distribution table. Notice that the degrees of freedom is \(n-1\).
  5. Step 5: If \(\chi^2 > \chi^2_{\alpha,n-1}\), reject the null hypothesis that the numbers are uniformly distributed. Otherwise, fail to reject the null hypothesis.

9.6 Auto-Correlation Test

This test is used to check the independence of the generated random numbers.

To study the autocorrelation, you first choose a lag \(m\) and you compute the autocorrelation for every \(m\) sample (\(R_i, R_{i+m}, R_{i+2m}, \ldots, R_{i+(M+1)m}\) where \((M+1)m < n\)).

The autocorrelation is computed as follows:

\[ \hat{\rho}_{im} = \frac{1}{M+1} \left[ \sum_{k=0}^M R_{i+km} R_{i+(k+1)m} \right] -0.25 \]

\[ \sigma_{\hat{\rho}_{im}} = \frac{\sqrt{13M+7}}{12(M+1)} \]

\[ Z_0 = \frac{\hat{\rho}_{im}}{\sigma_{\hat{\rho}_{im}}} \]

If \(-Z_{\alpha/2} \leq Z_0 \leq Z_{\alpha/2}\), then we fail to reject the null hypothesis that the numbers are independent. Otherwise, we reject the null hypothesis.