Monte Carlo Simulation Methodology

To estimate $\theta = E[g(X)]$ where $X \sim f$, the Monte Carlo estimator is:

$$\hat{\theta}_N = \frac{1}{N} \sum_{i=1}^{N} g(X_i), \quad X_i \sim f \text{ i.i.d.}$$

Properties:

• Unbiased: $E[\hat{\theta}_N] = \theta$
• Consistent: $\hat{\theta}_N \xrightarrow{a.s.} \theta$ (law of large numbers)
• CLT: $\sqrt{N}(\hat{\theta}_N - \theta) \xrightarrow{d} \mathcal{N}(0, \sigma_g^2)$

For BER estimation: $g(X_i) = \mathbf{1}[\text{bit } i \text{ is in error}]$ and $\hat{P}_b = k/N$ where $k$ is the error count.

The canonical BER simulation has this structure:

def simulate_ber(snr_db, n_bits, modulation='bpsk'):
    snr_lin = 10 ** (snr_db / 10)
    rng = np.random.default_rng(42)

    # 1. Generate random bits
    bits = rng.integers(0, 2, size=n_bits)

    # 2. Modulate: BPSK -> {-1, +1}
    x = 2 * bits - 1

    # 3. Channel: AWGN
    noise_std = 1 / np.sqrt(2 * snr_lin)
    y = x + noise_std * rng.standard_normal(n_bits)

    # 4. Detect: hard decision
    bits_hat = (y > 0).astype(int)

    # 5. Count errors
    n_errors = np.sum(bits != bits_hat)
    ber = n_errors / n_bits
    return ber, n_errors

Key rules:

• Process bits in large blocks (vectorized), not one at a time
• Use rng objects, not np.random.seed()
• Track error count separately for confidence intervals

Two common stopping criteria for BER simulations:

Fixed-count stopping: Run until $k_{\min}$ errors are observed:

n_errors = 0
n_bits_total = 0
while n_errors < k_min:
    ber_block, errors_block = simulate_block(snr_db, block_size)
    n_errors += errors_block
    n_bits_total += block_size
ber = n_errors / n_bits_total

Precision-based stopping: Run until the relative CI width $< \epsilon$:

while ci_width / ber_hat > epsilon:
    # ... run more blocks

Typical choices: $k_{\min} = 100$ errors, or $\epsilon = 10\%$ relative precision at 95% confidence.

The bit energy to noise ratio $E_b/N_0$ and the SNR per symbol are related by:

$$\frac{E_b}{N_0} = \frac{\mathrm{SNR}}{R \cdot \log_2 M}$$

where $R$ is the code rate and $M$ is the modulation order.

def ebno_to_snr(ebno_db, bits_per_symbol, code_rate=1.0):
    return ebno_db + 10 * np.log10(bits_per_symbol * code_rate)

Monte Carlo is embarrassingly parallel: each trial is independent. Use multiprocessing.Pool or concurrent.futures.ProcessPoolExecutor to distribute SNR points across CPU cores:

from concurrent.futures import ProcessPoolExecutor
from numpy.random import SeedSequence

def worker(args):
    snr_db, seed = args
    rng = np.random.default_rng(seed)
    return simulate_ber(snr_db, n_bits=1_000_000, rng=rng)

ss = SeedSequence(42)
seeds = ss.spawn(len(snr_range))
with ProcessPoolExecutor(max_workers=8) as pool:
    results = list(pool.map(worker, zip(snr_range, seeds)))