Simulating Exponential Random Variates using Julia

This demonstration provides a more comprehensive and practical example of how to simulate exponential random numbers in Julia. It covers both common methods, includes verification steps, and offers optional visualization to help you understand and trust the results. Remember to install the Plots package if you intend to use the histogram functionality.

Method 1: Using the built-in rand function with the Exponential distribution

# Method 1: Using the built-in `rand` function with the Exponential distribution

# Specify the rate parameter (lambda).  The mean of the Exponential distribution is 1/lambda.
lambda = 2.0  # Example rate

# Generate a single exponential random number
exponential_rn = rand(Exponential(1/lambda))
println("Single Exponential random number (Method 1): ", exponential_rn)

# Generate an array of exponential random numbers
num_samples = 1000
exponential_rns = rand(Exponential(1/lambda), num_samples)

# Verify the mean and variance (should be approximately 1/lambda and 1/lambda^2 respectively)
mean_simulated = mean(exponential_rns)
var_simulated = var(exponential_rns)
println("Mean of simulated (Method 1): ", mean_simulated)
println("Variance of simulated (Method 1): ", var_simulated)

Method 2: Using the inverse transform sampling method


# Method 2: Using the inverse transform sampling method

# Generate a uniform random number between 0 and 1
uniform_rn = rand()

# Apply the inverse of the CDF of the exponential distribution: -log(1-U)/lambda
exponential_rn_inverse = -log(1 - uniform_rn) / lambda
println("Single Exponential random number (Method 2): ", exponential_rn_inverse)


# Generate an array of exponential random numbers using inverse transform sampling
uniform_rns = rand(num_samples)
exponential_rns_inverse = -log.(1 .- uniform_rns) ./ lambda  # Use . for element-wise operations

# Verify the mean and variance
mean_simulated_inverse = mean(exponential_rns_inverse)
var_simulated_inverse = var(exponential_rns_inverse)
println("Mean of simulated (Method 2): ", mean_simulated_inverse)
println("Variance of simulated (Method 2): ", var_simulated_inverse)

Visualization the Distribution


# Plotting a histogram to visualize the distribution (optional, requires Plots.jl)
using Plots

histogram(exponential_rns, bins = 50,  xlabel = "Random Number", ylabel = "Frequency",
          title = "Histogram of Exponential Random Numbers (Method 1)",
          normalize = true, label = "Simulated") # Normalize to probability density

# Overlay the theoretical PDF for comparison
x = range(0, 5/lambda, length=100) # Adjust range as needed
plot!(x, lambda * exp.(-lambda*x), label = "Theoretical PDF", linewidth=2)

savefig("exponential_histogram_method1.png")


histogram(exponential_rns_inverse, bins = 50,  xlabel = "Random Number", ylabel = "Frequency",
          title = "Histogram of Exponential Random Numbers (Method 2)",
          normalize = true, label = "Simulated") # Normalize to probability density

# Overlay the theoretical PDF for comparison
x = range(0, 5/lambda, length=100) # Adjust range as needed
plot!(x, lambda * exp.(-lambda*x), label = "Theoretical PDF", linewidth=2)

savefig("exponential_histogram_method2.png")

Explanation and Key Improvements:

  1. Two Methods: The code demonstrates two common ways to generate exponential random numbers:

    • Method 1: Using the built-in rand(Exponential(1/lambda)) function. This is the most straightforward and efficient way in Julia. Note that Exponential(θ) in Julia uses the mean (θ) parameterization, which is the inverse of the rate (λ) often used in other contexts.
    • Method 2: Using the inverse transform sampling method. This method is more general and can be applied to other distributions as well. It involves generating a uniform random number and then transforming it using the inverse of the cumulative distribution function (CDF) of the exponential distribution.

  2. Clearer Parameterization: The code emphasizes the use of the rate parameter (lambda) and clearly explains the relationship between the rate and the mean of the exponential distribution (mean = 1/lambda).

  3. Verification: The code calculates and prints the mean and variance of the generated random numbers to verify that they are close to the expected theoretical values (1/lambda and 1/lambda^2, respectively). This is a good practice to ensure the simulation is working correctly.

  4. Histograms (Optional): The code includes an optional section that uses the Plots.jl package to create histograms of the generated random numbers and overlay the theoretical probability density function (PDF). This allows you to visually compare the simulated distribution to the expected distribution. You’ll need to install the Plots package if you want to use this part: ] add Plots in the Julia REPL. The savefig function is used to save the plot as a PNG image.

  5. Element-wise operations: The code uses the dot operator (.) for element-wise operations when generating arrays of random numbers using the inverse transform method. This is crucial for efficiency in Julia.

  6. Comments and Clarity: The code is well-commented to explain each step, making it easier to understand and modify.