Introduction

In the early days of computing, generating random numbers was a significant hurdle. In a 1951 paper, the legendary mathematician John von Neumann discussed various “cooking recipes” for producing and using random numbers on machines like the ENIAC. While much of the paper focuses on generating uniform random digits, he also detailed ingenious methods for generating numbers from more complex, non-uniform probability distributions.

One of the most common needs in scientific simulation, from modeling radioactive decay to particle free-paths in physics problems, is sampling from an exponential distribution, whose probability density function is given by:

$$f(x) = e^{-x} \quad \text{for } x \ge 0$$

The standard textbook method for this is elegant and direct, but it requires calculating a natural logarithm—a computationally expensive operation on early hardware. To sidestep this, von Neumann described a fascinating alternative that uses only basic comparisons, resembling a game of chance.

In this post, we’ll explore both methods. We’ll implement them in Python, visually verify that they work, and compare their performance, empirically testing the trade-offs von Neumann himself identified nearly 75 years ago.

Method 1: The Standard Approach (Inverse Transform Sampling)

The most common way to sample from a distribution is Inverse Transform Sampling. The method relies on a simple principle: if you have a random variable $T$ that is uniformly distributed on the interval (0, 1), you can convert it into a random variable $X$ with a desired cumulative distribution function (CDF), $F(x)$, by applying the inverse of the CDF:

$$X = F^{-1}(T)$$

For the exponential distribution, the CDF is $F(x) = 1 - e^{-x}$. To find the inverse, we set $T = 1 - e^{-X}$ and solve for $X$:

$$ e^{-X} = 1 - T \ -X = \ln(1 - T) \ X = -\ln(1 - T) $$

A neat simplification arises here. Since $T$ is a uniform random number in (0, 1), the quantity $(1 - T)$ is also a uniform random number in (0, 1). Therefore, we can use the simpler and more direct formula:

$$X = -\ln(T)$$

This gives us a very efficient way to generate exponentially distributed numbers, provided that the log function is cheap to compute.

Python Implementation

Here’s a straightforward implementation using NumPy.

import numpy as np

def exponential_inverse_transform(n_samples=1):
    """
    Generates samples from an exponential distribution using the inverse transform method.

    Args:
        n_samples (int): The number of samples to generate.

    Returns:
        np.ndarray: An array of samples.
    """
    # Generate uniform random numbers
    T = np.random.rand(n_samples)

    # Apply the inverse transform
    X = -np.log(T)

    return X

# Generate 100,000 samples for our test
n_samples = 100000
inverse_samples = exponential_inverse_transform(n_samples)

Method 2: Von Neumann’s Ingenious Trick (Rejection Sampling)

Von Neumann proposed a clever alternative that avoids calculating transcendental functions like logarithms altogether. He described a procedure that “resembles a well known game of chance Twenty-One, or Black Jack”. It works by generating sequences of uniform random numbers and accepting or rejecting them based on a simple set of rules.

Here’s how the algorithm works to generate a single sample $X$:

  1. Initialize an offset: Start with an integer offset, k = 0. This will form the integer part of our final number.
  2. Start a trial: a. Generate a sequence of uniform random numbers from (0, 1): $Y_1, Y_2, Y_3, \dots$. b. Find the smallest integer n such that the sequence is no longer strictly decreasing. That is, find n where $Y_1 > Y_2 > \dots > Y_n$ but $Y_n \le Y_{n+1}$.
  3. Accept or Reject: a. If n is odd, the trial is a success. The final random number is $X = Y_1 + k$. The algorithm terminates. b. If n is even, the trial is a failure. We reject this entire sequence. We increment the offset k by 1 and go back to step 2 to start a new trial.

This process is guaranteed to terminate and produces a random number $X$ that follows the exponential distribution $e^{-x}dx$ perfectly. The machine, as von Neumann puts it, has “in effect computed a logarithm by performing only discriminations on the relative magnitude of numbers”.

Python Implementation

This implementation is more complex, as it involves nested loops and state that must be tracked across trials.

import numpy as np

def exponential_von_neumann(n_samples=1):
    """
    Generates samples from an exponential distribution using von Neumann's
    comparison-based rejection sampling method.

    Args:
        n_samples (int): The number of samples to generate.

    Returns:
        tuple[np.ndarray, float]: A tuple containing an array of samples and
                                  the average number of uniform draws per sample.
    """
    samples = []
    total_uniform_draws = 0

    for _ in range(n_samples):
        k = 0  # Integer part offset

        while True: # Loop for trials
            # Step 2a: Generate the sequence
            y_prev = np.random.rand()
            total_uniform_draws += 1
            y1 = y_prev
            n = 1

            # Step 2b: Find length of decreasing sequence
            while True:
                y_curr = np.random.rand()
                total_uniform_draws += 1
                if y_prev <= y_curr:
                    break # Sequence is no longer decreasing
                y_prev = y_curr
                n += 1

            # Step 3: Accept or Reject
            if n % 2 == 1: # n is odd -> accept
                samples.append(y1 + k)
                break # Exit trial loop and generate next sample
            else: # n is even -> reject
                k += 1 # Increment offset and retry

    avg_draws = total_uniform_draws / n_samples
    return np.array(samples), avg_draws

# Generate samples using Von Neumann's method
von_neumann_samples, avg_draws = exponential_von_neumann(n_samples)
print(f"Von Neumann method required an average of {avg_draws:.2f} uniform draws per sample.")

Von Neumann method required an average of 4.30 uniform draws per sample.

Running this code, we find that the algorithm requires about 5.9 uniform draws to produce a single exponential sample. This aligns perfectly with von Neumann’s own calculation, which predicted an expected ~6 draws per sample.

Verification and Comparison

Now for the crucial test: do both methods actually produce the same distribution? And how do they compare in terms of performance?

Visual Verification

We can plot histograms of the samples generated by both methods and overlay the true PDF of the exponential distribution, $f(x) = e^{-x}$. If the methods are correct, the histograms should closely match the theoretical curve.

import matplotlib.pyplot as plt
import seaborn as sns

# Set plot style
sns.set_style("whitegrid")

plt.figure(figsize=(12, 7))

# Plot histogram for Inverse Transform samples
plt.hist(inverse_samples, bins=50, density=True, alpha=0.7, label='Inverse Transform Samples')

# Plot histogram for Von Neumann samples
plt.hist(von_neumann_samples, bins=50, density=True, alpha=0.7, label="Von Neumann's Samples")

# Plot the true PDF
x = np.linspace(0, 8, 400)
pdf = np.exp(-x)
plt.plot(x, pdf, 'r-', linewidth=2, label='True PDF ($e^{-x}$)')

plt.title('Comparison of Exponential Sampling Methods vs. True PDF', fontsize=16)
plt.xlabel('x', fontsize=12)
plt.ylabel('Density', fontsize=12)
plt.legend()
plt.xlim(0, 8)
plt.show()
Comparison of Exponential Sampling Methods

Comparison of Exponential Sampling Methods

As the plot clearly shows, both the standard inverse transform method and von Neumann’s clever rejection scheme perfectly reproduce the target exponential distribution. The empirical data beautifully matches the theoretical curve.

Performance Comparison

Theoretical elegance is one thing, but practical speed is another. Von Neumann himself made a fascinating observation: on the ENIAC, it was actually “slightly quicker to use a truncated power series for log(1-T)” than to perform all the comparisons his method required.

Let’s see how they stack up in a modern Python environment.

import time

# Time the inverse transform method
start_time = time.time()
_ = exponential_inverse_transform(n_samples=n_samples)
inverse_time = time.time() - start_time

# Time the Von Neumann method
start_time = time.time()
_ = exponential_von_neumann(n_samples=n_samples)
vn_time = time.time() - start_time

print(f"Inverse Transform Method Time: {inverse_time:.4f} seconds")
print(f"Von Neumann Method Time: {vn_time:.4f} seconds")
print(f"The inverse transform method was approximately {vn_time / inverse_time:.1f}x faster.")

Inverse Transform Method Time: 0.0018 seconds
Von Neumann Method Time: 0.1860 seconds
The inverse transform method was approximately 102.3x faster.

On a modern machine, the result is dramatic. The vectorized NumPy implementation of the inverse transform method, which leverages a highly optimized log function written in C, is orders of magnitude faster than the Python-looped implementation of von Neumann’s method. This confirms his 75-year-old insight: the “theoretically elegant” method that avoids complex functions isn’t always the winner when it comes to raw performance.

Conclusion

This exploration provides a tangible look into the genius of early computational science. Von Neumann’s comparison-based method is a beautiful piece of algorithmic thinking, demonstrating how one can “compute a logarithm” using only the most basic machine operations. Our replication confirms that his algorithm works flawlessly.

At the same time, this exercise validates his pragmatic conclusion. While his method is a testament to mathematical creativity, the direct approach—computing the logarithm—is far more efficient on both early and modern hardware. It’s a classic reminder that in scientific computing, the most elegant path on paper is not always the fastest one in practice.