Exercise

Generate a random variable generator using the inversion method for the following probability density function:

\[f(x) = 2.5x \cdot \sqrt{x} = 2.5x^{3/2} \quad \text{for } 0 < x < 1\]

Otherwise \(f(x) = 0\).

Step 1: Verify that f(x) is a probability density function

A function \(f(x)\) is a valid PDF if:

\[\int_{-\infty}^{+\infty} f(x)dx = 1\]

In our case:

\[\int_{0}^{1} 2.5x^{3/2} dx = 2.5 \cdot \left[\frac{2}{5}x^{5/2}\right]_0^1 = 2.5 \cdot \frac{2}{5} = 1 \quad \]

Step 2: Calculate F(x)

The CDF is defined as:

\[F(x) = \int_{0}^{x} f(t)dt = \int_{0}^{x} 2.5t^{3/2} dt\]

Solving the integral:

\[F(x) = 2.5 \cdot \left[\frac{2}{5}t^{5/2}\right]_0^x = 2.5 \cdot \frac{2}{5} \cdot x^{5/2} = x^{5/2}\]

Therefore: \(F(x) = x^{5/2}\) for \(0 \leq x \leq 1\).

Step 3: Find F^(-1)(u)

To apply the inversion method, we need to find the inverse of F(x).

Given \(F(x) = x^{5/2}\), we solve for x:

\[u = x^{5/2}\]

\[x = u^{1/(5/2)} = u^{2/5}\] Therefore: \(F^{-1}(u) = u^{2/5}\) for \(0 \leq u \leq 1\).

Visualization

# Define the CDF
F_x <- function(x) {
  x^(5/2)
}

# Define the inverse CDF
F_inv <- function(u) {
  u^(2/5)
}

# Create sequences
x_vals <- seq(0, 1, length.out = 100)
F_vals <- F_x(x_vals)

u_vals <- seq(0, 1, length.out = 100)
x_inv <- F_inv(u_vals)

# Plot both F(x) and F^(-1)(u)
par(mfrow = c(1, 2))

# Plot F(x)
plot(x_vals, F_vals, 
     type = "l", 
     col = "blue", 
     lwd = 2,
     main = "F(x) = x^(5/2)",
     xlab = "x",
     ylab = "F(x)",
     las = 1)
grid()

# Plot F^(-1)(u)
plot(u_vals, x_inv, 
     type = "l", 
     col = "red", 
     lwd = 2,
     main = "F⁻¹(u) = u^(2/5)",
     xlab = "u",
     ylab = "x = F⁻¹(u)",
     las = 1)
grid()

Step 3: Generator Implementation

Now we implement the generator using the inversion method. The algorithm is:

  1. Generate \(U \sim \text{Uniform}(0,1)\)
  2. Compute \(X = F^{-1}(U) = U^{2/5}\)
  3. \(X\) follows the distribution with density \(f(x) = 2.5x^{3/2}\)
# Set seed for reproducibility
set.seed(27)

# Number of simulations
nsim <- 10000

# Generator function using inversion method
simRV <- function() {
  u <- runif(1)
  x <- u^(2/5)
  return(x)
}

# Generate nsim random values
RV_simulated <- replicate(nsim, simRV())

# First 10 values
cat("First 10 simulated values:\n")
## First 10 simulated values:
print(head(RV_simulated, 10))
##  [1] 0.98860286 0.37085958 0.94749890 0.64121029 0.54796983 0.69428589
##  [7] 0.35005127 0.09029901 0.45165558 0.51669968
# Summary
summary(RV_simulated)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01305 0.57105 0.75611 0.71163 0.88901 0.99997

Step 4: Validation

We validate our generator by comparing the simulated histogram with the theoretical PDF.

# Histogram with theoretical PDF overlay
hist(RV_simulated, 
     breaks = 30, 
     freq = FALSE,
     col = "lightblue",
     border = "white",
     main = "Simulated Distribution vs Theoretical PDF",
     xlab = "x",
     ylab = "Density",
     las = 1,
     xlim = c(0, 1))

# Add theoretical PDF curve
curve(2.5 * x^(3/2), 
      from = 0, 
      to = 1, 
      add = TRUE, 
      col = "red", 
      lwd = 2)

# Add legend
legend("topleft", 
       legend = c("Simulated", "Theoretical f(x)"),
       fill = c("lightblue", NA),
       border = c("black", NA),
       lty = c(NA, 1),
       col = c(NA, "red"),
       lwd = c(NA, 2),
       bty = "n")

# Summary statistics
cat("Theoretical mean: E[X] = ", integrate(function(x) x * 2.5 * x^(3/2), 0, 1)$value, "\n")
## Theoretical mean: E[X] =  0.7142857
cat("Simulated mean:", mean(RV_simulated), "\n")
## Simulated mean: 0.7116339
cat("Theoretical variance:", 
    integrate(function(x) x^2 * 2.5 * x^(3/2), 0, 1)$value - 
    (integrate(function(x) x * 2.5 * x^(3/2), 0, 1)$value)^2, "\n")
## Theoretical variance: 0.04535147
cat("Simulated variance:", var(RV_simulated), "\n")
## Simulated variance: 0.04581319

Conclusions

The inversion method has been successfully implemented to generate random variables from the density function \(f(x) = 2.5x^{3/2}\).

The graphical comparison shows excellent agreement between the simulated histogram and the theoretical PDF curve. The simulated mean and variance are very close to their theoretical values, confirming the correctness of our implementation.