An Original Description of the Problem

Generate n = 10,000 values for a random variable X that has the probability density function f(x) = (x^2)/9, for 0 < x < 3, using the rejection method. Create a density histogram for the generated values, superimpose the probability density function of the target distribution on the histogram, estimate the mean, standard deviation, and the probability that X < 2 using the generated values, and compare these estimated values with the theoretical ones. Report the rejection rate.

Simulated Data

Data from print(head(x_values)): 2.8214019, 2.8705000, 0.7382632, 2.6686179, 2.8890727, 2.2753786

Summary Graphs

Findings Discussion

The simulation using the rejection method to generate random variables X with the probability density function f(x) = (x^2)/9 over 0 < x < 3 shows high accuracy. The estimated mean (2.2546) and standard deviation (0.5767) are very close to the theoretical values (2.25 and 0.5809, respectively). The estimated probability P(X < 2) is 0.2936, matching well with the theoretical probability of 0.2963. The histogram with the target PDF superimposed visually confirms this accuracy. Although the rejection rate is high at about 66.92%, the method is simple and effective for generating samples from complex distributions. The small differences between estimated and theoretical values highlight the method’s reliability.

Code

n <- 10000 f <- function(x) (x^2) / 9 g <- function(x) 1 / 3 M <- 3

set.seed(123) x_values <- numeric(n) i <- 1 rejections <- 0

while(i <= n) { u <- runif(1, 0, 3) v <- runif(1, 0, 1) if (v <= f(u) / (M * g(u))) { x_values[i] <- u i <- i + 1 } else { rejections <- rejections + 1 } }

hist(x_values, probability = TRUE, breaks = 30, main = “Density Histogram of Generated Values”, xlab = “x”, ylab = “Density”)

curve((x^2) / 9, add = TRUE, col = “red”, lwd = 2)

legend(“topleft”, legend = c(“Target PDF”), col = c(“red”), lty = 1, lwd = 2)

estimated_mean <- mean(x_values) estimated_sd <- sd(x_values) estimated_prob <- mean(x_values < 2)

theoretical_mean <- integrate(function(x) x * (x^2) / 9, 0, 3)\(value theoretical_sd <- sqrt(integrate(function(x) (x^2) * (x^2) / 9, 0, 3)\)value - theoretical_mean^2) theoretical_prob <- integrate(function(x) (x^2) / 9, 0, 2)$value

print(paste(“Estimated Mean:”, estimated_mean)) print(paste(“Theoretical Mean:”, theoretical_mean)) print(paste(“Estimated Standard Deviation:”, estimated_sd)) print(paste(“Theoretical Standard Deviation:”, theoretical_sd)) print(paste(“Estimated P(X < 2):”, estimated_prob)) print(paste(“Theoretical P(X < 2):”, theoretical_prob))

rejection_rate <- rejections / (rejections + n) print(paste(“Rejection Rate:”, rejection_rate))

print(head(x_values))