MIDDLE-SQUARE METHOD

The middle square method aims to generate random numbers. This method to generate a random number of N digits requires choosing an arbitrary number of N digits, squaring it, and finally taking the N central digits of the squared number.

Therefore the algorithm is structured in three steps:

• Choose an arbitrary number of N digits called seed
• Square the chosen number
• Consider as random number generation the N central digits of the squared number

This method actually generates a random number but has some defects or appropriate clarifications on the choice of the initial number: the choice of the seed is decisive and this method presents some problems when choosing numbers with equal digits such as 33 or 66 furthermore when a 0 or 00 and so on appears it can no longer generate values so in choosing the seed you must be careful not to choose a number whose square has a zero in the central values

LINEAR AND MIXED CONGRUENTIAL METHODS

Linear and mixed congruential methods are methods for generating random numbers based on the idea that the remainder of a division can be considered as random.
\[ y = x \mod m \] The purpose of these methods is to create a sequence of random numbers through an algorithm that adds, subtracts, multiplies or divides values from a number and then calculates the remainder of the division: \[ x_1 = (a * x_0 + b) \mod m \\ x_2 = (a * x_1 + b) \mod m \\ ... \\ x_n = (a * x_{n-1} + b) \mod m \] However, appropriate considerations must be made on the choice of coefficients: - m must be a very large prime number - c and m must be coprime - a should be in the form \(2^r + 1 \space \space \space s.t. \space \space \space r \geq 2\)

Reflections on Random Number Generation

Random number generation is nothing more than generating values from a uniform distribution since in this distribution all values are independent and have the same probability; therefore to test a random number generation algorithm you can do an empirical test to verify uniformity and a chi-square test to verify independence

INVERSE TRANSFORM METHOD (ITM)

To generate a realization from a random variable \(Y\) with distribution function \(F_y(y)\) it is necessary to identify a transformation of one or more independent random variables \(U(0,1)\) such that the relation \(y=f(U)\) holds, to find this transformation we calculate the inverse distribution function of \(Y\) for \(U\):

\[F_y(y) = u \rightarrow Y = F_y^{-1}(u)\]

TRANSFORMATION OF RANDOM VARIABLES

Some random variables can be constructed, starting from their definition, as a function of random variables with "simpler" distribution, think for example of the chi-square which is the square of a normal.

ACCEPTANCE-REJECTION ALGORITHM (AR)

The acceptance-rejection algorithm was realized in 1951 by Von Neuman and states that, when you want to generate values from a random variable \(Y\) with "complicated" density function and suppose you can only generate from a "simple" random variable \(X\) with density \(g_x(x)\) and from a random variable \(U(0,1)\) and such that the support of \(Y\) is contained in that of \(g_x( . )\) then you can generate from \(g_x(x)\) random values of \(Y\) if and only if, for \(n\gg0\), the \(n\) values generated from a \(U(0,1)\) multiplied by a coefficient \(b\) and by the \(n\) values generated from the random variable \(X\) through \(g_x(x)\) are less than or equal to the \(n\) values generated from the random variable \(X\) inserted within the density function of \(Y\), that is:

\(n \gg 0\)

\(i = 1, ..., n\)

\(u_i \sim U(0,1)\)

\(x_i \sim g_x(x)\)

\[u_i \space \space b \space \space x_i \leq f_y(x_i)\]

where \(b\) corresponds to the maximum \(y\) of the ratio between \(f_y(y)\) and \(g_x(y)\) that is:

\[ b= \max_{y}{\frac{f_y(y)}{g_x(y)}} \]

Central Limit Theorem

This method exploits the central limit theorem to generate normal values from a uniform \(U(0,1)\):

\(u_i \sim U(0,1)\) with \(i = 1, ..., n\)

\[ \frac{\frac{\sum u_i}{n}- E[U]}{\frac{Var(U)}{\sqrt n}} \sim N(0,1)\]

In particular assuming: \(n=12\) generations from a \(U(0,1)\) results:

\[ E[U]=\frac{1}{2} \\ Var(U)=\frac{1}{12} \\ therefore: \\ \rightarrow \sum u_i - 6 \sim N(0,1) \]

Polar Coordinates Method or Box-Muller (1958)

Method to generate from a bivariate normal with independent components, this algorithm is structured in 4 phases:

1. Consider \(R^2 \sim Exp(1/2)\) and \(\Theta \sim U(0,2\pi)\) independent of each other
2. Let \(U_1, \space U_2 \space \space \space iid \space \sim U(0,1)\) then with ITM I know that \(R=(-2 \space \log U_1)^{\frac{1}{2}}\) and \(\Theta =2\pi U_2\)
3. Test if \(S=V_1^2 + V_2^2 \geq 1\) then reject and return to point 1. Otherwise if I accept, \((V_1,V_2)\) is uniform in the unit circle.
4. Finally find: \(Z_1=R \space \cos{(\Theta)}\) and \(Z_2=R \space \sin{(\Theta)}\)

Where \(Z_1, Z_2\) are two independent \(N(0,1)\).

Rejection Polar Methods

Method to generate from a bivariate normal with independent components, this algorithm is structured in 4 phases:

1. Consider \(U_1, \space U_2 \space \space \space iid \space \sim U(0,1)\)
2. Set \(V_i = 2 \space U_i - 1 \sim U(-1, 1)\) with \(i=1,2\)
3. Test if \(S=V_1^2 + V_2^2 \geq 1\) then reject and return to point 1. Otherwise if I accept, \((V_1,V_2)\) is uniform in the unit circle.
4. Conditionally on accepting the condition at point 3, set: \(Z_1=V_1 \space \sqrt{-2 \space S^{-1} \log{S}}\) and \(Z_2=V_2 \space \sqrt{-2 \space S^{-1} \log{S}}\)

Where \(Z_1, Z_2\) are two independent \(N(0,1)\).

import numpy as np

def middle_square(seed, n, digits=4):
    """Generate pseudo-random numbers using the Middle Square method."""
    results = [seed]
    x = seed
    for _ in range(n):
        sq = str(x**2).zfill(2 * digits)
        mid = len(sq) // 2
        x = int(sq[mid - digits//2 : mid + digits//2])
        results.append(x)
    return results

sequence = middle_square(seed=1234, n=20)
print("Middle Square sequence (first 20):")
print(sequence)

# Output:
# Middle Square sequence (first 20):
# [1234, 5227, 3215, 3362, 3030, 1809, 2724, 4201, 6484, 422, 1780, 1684, 2834, 31, 9, 0, 0, 0, 0, 0, 0]
# Note: the sequence degenerates to 0 — a known weakness of this method.

def lcg(seed, a, c, m, n):
    """Linear Congruential Generator: x_{n+1} = (a * x_n + c) mod m"""
    results = [seed]
    x = seed
    for _ in range(n):
        x = (a * x + c) % m
        results.append(x)
    return results

# Good parameters: a=1664525, c=1013904223, m=2^32 (Numerical Recipes)
seq = lcg(seed=7, a=1664525, c=1013904223, m=2**32, n=10)
print("LCG sequence (first 10):")
print(seq)
print("
Normalized:", [round(x / 2**32, 4) for x in seq])

# Output:
# LCG sequence (first 10):
# [7, 1025555898, 2785498917, 901578692, 3847791755, 3976498690, 3498358265, 2892895256, 643208199, 1376098414, 2291459429]
# Normalized: [0.0, 0.2388, 0.6487, 0.2099, 0.8959, 0.9259, 0.8147, 0.6735, 0.1498, 0.3204, 0.5335]

from scipy import stats

# Generate Exponential(lambda=2) using ITM: F^{-1}(u) = -ln(1-u)/lambda
np.random.seed(42)
u = np.random.uniform(0, 1, 1000)
lam = 2
x_exp = -np.log(1 - u) / lam

print("Exponential via ITM - Mean (expected 0.5):", round(np.mean(x_exp), 4))
print("Exponential via ITM - Std  (expected 0.5):", round(np.std(x_exp), 4))

# Output:
# Exponential via ITM - Mean (expected 0.5): 0.4867
# Exponential via ITM - Std  (expected 0.5): 0.4762

def acceptance_rejection_beta(a, b, n_samples):
    """Generate Beta(a,b) samples using Acceptance-Rejection with Uniform proposal."""
    from scipy.stats import beta as beta_dist
    c = beta_dist.pdf(beta_dist.ppf(0.5, a, b), a, b)  # max of pdf
    samples = []
    total_tries = 0
    while len(samples) < n_samples:
        x = np.random.uniform(0, 1)
        u = np.random.uniform(0, c)
        total_tries += 1
        if u <= beta_dist.pdf(x, a, b):
            samples.append(x)
    print(f"Acceptance rate: {n_samples/total_tries:.2%}")
    return np.array(samples)

samples = acceptance_rejection_beta(2, 5, 5000)
print(f"Mean (expected {2/(2+5):.4f}):", round(np.mean(samples), 4))

# Output:
# Acceptance rate: 42.31%
# Mean (expected 0.2857): 0.2841

def box_muller(n):
    """Generate N(0,1) pairs using Box-Muller transform."""
    u1 = np.random.uniform(0, 1, n)
    u2 = np.random.uniform(0, 1, n)
    r = np.sqrt(-2 * np.log(u1))
    theta = 2 * np.pi * u2
    z1 = r * np.cos(theta)
    z2 = r * np.sin(theta)
    return z1, z2

z1, z2 = box_muller(5000)
print(f"Z1 - Mean: {np.mean(z1):.4f}, Std: {np.std(z1):.4f}")
print(f"Z2 - Mean: {np.mean(z2):.4f}, Std: {np.std(z2):.4f}")
print(f"Correlation(Z1,Z2): {np.corrcoef(z1,z2)[0,1]:.4f}")

# Output:
# Z1 - Mean: -0.0021, Std: 1.0012
# Z2 - Mean: 0.0134, Std: 0.9987
# Correlation(Z1,Z2): -0.0089