Website link : https://rpubs.com/Isaiah-Mireles/1446105

Q1)

(a) Write R code using the inverse transform method to generate random numbers from the distribution of X.

# random samples
n <- 10000

# Generate Uniform(0,1) random numbers
u <- runif(n)

# Allocate memory
x <- numeric(n)

# Inverse transform method
x[u <= 0.30] <- 1
x[u > 0.30 & u <= 0.50] <- 2
x[u > 0.50 & u <= 0.65] <- 3
x[u > 0.65 & u <= 0.75] <- 4
x[u > 0.75] <- 5

table(x) / n

## x
##      1      2      3      4      5 
## 0.2971 0.1937 0.1576 0.0971 0.2545

(b) Generate 10,000 random numbers and draw a bar chart.

(table(x) / n) |> barplot()

(c) Compare the sample relative frequencies with the theoretical probability distribution.

true <- c(0.3, 0.2, 0.15, 0.1, 0.25)

(table(x) / n) - true

## x
##       1       2       3       4       5 
## -0.0029 -0.0063  0.0076 -0.0029  0.0045

as we can see the prob. are similar, very minimal difference – just sampling error

Q2) Suppose we want to generate random samples from a Poisson distribution with parameter \(\lambda\) using the inverse CDF method.

a) Design an algorithm to generate random numbers from the Poisson distribution with parameter \(\lambda\) using the inverse CDF method.

Below I describe the recursive algorithm for poiss() :

b-c) Apply Algorithm – here we apply recursive func.

set.seed(123)
n <- 10000
lambda <- 4.2

x <- numeric(n)

for (i in 1:n) {
  u <- runif(1) # generate a single value

  xi <- 0
  p <- exp(-lambda) # x=0
  cdf <- p
  
  # stop when we reach the max value 
  while (cdf < u) {
    xi <- xi + 1
    p <- p * lambda / xi # recursive def.
    cdf <- cdf + p
  }

  x[i] <- xi
}

Compare :

hist(x,
     probability = TRUE,
     breaks = seq(-0.5, max(x) + 0.5, by = 1),
     col = "lightgray",
     main = "Inverse CDF Poisson Samples",
     xlab = "x")

k <- 0:max(x)
points(k, dpois(k, lambda), col = "red", pch = 19)
lines(k, dpois(k, lambda), col = "red", lwd = 2)

as we can see, we did a good job

Q3) Consider the probability mass function provided in Problem 1. Suppose we want to generate random samples from this distribution using the acceptance-rejection method.

(a)

lets suppose we just assume for each, the prob is equal 0.2

# Target distribution
x <- 1:5
p_x <- c(0.30, 0.20, 0.15, 0.10, 0.25)

# Proposal distribution
g <- rep(0.2, 5)

(b) Acceptance Rate

# Compute M
M <- max(p_x / g)

# Acceptance probabilities
p_accept <- p_x / (M * g)

# Number of proposals
n <- 10000

# Generate proposals from g : here we assign equal prob. g (0.2 each)
x_prop <- sample(x, size = n, replace = TRUE, prob = g)

# Independent uniforms
u <- runif(n)

# Accept/reject
accepted <- u <= p_accept[x_prop]

print("Empirical Acceptance Rate")

## [1] "Empirical Acceptance Rate"

mean(accepted)

## [1] 0.6667

print("Theoretical Acceptance Rate")

## [1] "Theoretical Acceptance Rate"

1/M

## [1] 0.6666667

pretty high and close to theoretical

# Accepted samples
x_accept <- x_prop[accepted]

print("est.")

## [1] "est."

table(x_accept) / length(x_accept)

## x_accept
##         1         2         3         4         5 
## 0.2977351 0.2083396 0.1486426 0.1022949 0.2429879

print("True prob.")

## [1] "True prob."

p_x

## [1] 0.30 0.20 0.15 0.10 0.25

not bad!
We may find better luck pushing certain x values more likely to match the true dist better

g <- c(.25, .2, .15, .1, .3)

# Compute M
M <- max(p_x / g)

# Acceptance probabilities
p_accept <- p_x / (M * g)

# Number of proposals
n <- 10000

# Generate proposals from g : 
x_prop <- sample(x, size = n, replace = TRUE, prob = g)

# Independent uniforms
u <- runif(n)

# Accept/reject
accepted <- u <= p_accept[x_prop]

print("Empirical Acceptance Rate")

## [1] "Empirical Acceptance Rate"

mean(accepted)

## [1] 0.8367

print("Theoretical Acceptance Rate")

## [1] "Theoretical Acceptance Rate"

1/M

## [1] 0.8333333

as you see, we did better

Suppose we identically match true

g <- p_x

# Compute M
M <- max(p_x / g)

# Acceptance probabilities
p_accept <- p_x / (M * g)

# Number of proposals
n <- 10000

# Generate proposals from g : 
x_prop <- sample(x, size = n, replace = TRUE, prob = g)

# Independent uniforms
u <- runif(n)

# Accept/reject
accepted <- u <= p_accept[x_prop]

print("Empirical Acceptance Rate")

## [1] "Empirical Acceptance Rate"

mean(accepted)

## [1] 1

print("Theoretical Acceptance Rate")

## [1] "Theoretical Acceptance Rate"

1/M

## [1] 1

we notice at this point, there is perfect performance on both

Here’s a more concise version:A proposal distribution (g(x)) that closely matches the target distribution (p(x)) is more efficient because it produces a higher acceptance rate and fewer rejected samples. In the ideal case where (g(x)=p(x)), every proposal is accepted, so both the theoretical and empirical acceptance rates are 100%. Using a uniform proposal is generally less efficient because it ignores the shape of (p(x)), leading to more rejections. The empirical acceptance rate may differ from the theoretical rate due to random sampling, but it converges to the theoretical value as the sample size increases.

Personal Curiosity

set.seed(123)

# Generate x values
x <- seq(0, 2*pi, length.out = 100)

# Random error ~ N(0, 1)
error <- rnorm(length(x), mean = 0, sd = 1)

# Example response: sin(x) + error
y <- sin(x) + error

# Combine into a data frame
df <- data.frame(x, y, error)

plot(x,y)
lines(x, sin(x), col = "red", lwd = 3)

# Fit models
m1 <- lm(y ~ x, data = df)
m2 <- lm(y ~ sin(x), data = df)

# Scatterplot
plot(df$x, df$y,
     pch = 1,
     xlab = "x",
     ylab = "y",
     main = "Comparison of Models")

# True model
lines(df$x, sin(df$x),
      col = "red",
      lwd = 3)

# m1: Linear regression
abline(m1,
       col = "blue",
       lwd = 2)

# m2: Sine regression
lines(df$x,
      predict(m2),
      col = "darkgreen",
      lwd = 2)

legend("topright",
       legend = c("True Model", "m1: y ~ x", "m2: y ~ sin(x)"),
       col = c("red", "blue", "darkgreen"),
       lwd = c(3, 2, 2),
       bty = "n", 
       cex = .5)