Website (suggested) version : https://rpubs.com/Isaiah-Mireles/1447001

1 Q1

A \(t\)-dist using the stochastic representation :

\[ T=\frac{Z}{\sqrt{V/\nu}}, \]

where \(Z\sim N(0,1)\), \(V\sim\chi^2_\nu\), and \(Z\) and \(V\) are independent.

Using the Box–Muller transform, let

\[ U_1,U_2 \stackrel{\text{iid}}{\sim} \operatorname{Unif}(0,1). \]

Then

\[ Z_1=\sqrt{-2\ln(U_1)}\cos(2\pi U_2), \]

and

\[ Z_2=\sqrt{-2\ln(U_1)}\sin(2\pi U_2) \]

are independent \(N(0,1)\) random variables. We repeatedly apply this transformation to generate \(11\) independent standard normal random variables, compute

\[ V=\sum_{i=2}^{11} Z_i^2 \sim \chi^2_{10}, \]

and obtain \(t(10)\) random samples using

\[ T=\frac{Z_1}{\sqrt{V/10}}. \]

set.seed(123)

n <- 10000
df <- 10

t_samples <- rep(NA, n)

for (sim in 1:n) {

  # Generate 11 independent standard normals
  normals <- rep(NA, df + 1)

  for (i in seq(1, df + 1, by = 2)) {

    U1 <- runif(1)
    U2 <- runif(1)

    R <- sqrt(-2 * log(U1))
    theta <- 2 * pi * U2

    normals[i] <- R * cos(theta)

    if (i + 1 <= df + 1) {
      normals[i + 1] <- R * sin(theta)
    }
  }

  # First normal is the numerator
  Z <- normals[1]

  # Remaining 10 normals form a chi-square(10)
  V <- sum(normals[2:(df + 1)]^2)

  # t-distribution sample
  t_samples[sim] <- Z / sqrt(V / df)
}

# Plt
hist(t_samples,
     probability = TRUE,
     breaks = 50,
     col = "lightblue",
     border = "black",
     main = "10,000 Samples from t(10)",
     xlab = "t",
     ylab = "Density")

# Overlay the true t(10) density
curve(dt(x, df = df),
      add = TRUE,
      col = "red",
      lwd = 2)

2 Q2

2.1 (a)

Let

\[ X \sim \chi^2_1, \]

and define

\[ Y = \sqrt{X}. \]

Since \(X = Y^2\), the change-of-variables formula gives

\[ f_Y(y) = f_X(y^2) \left| \frac{d}{dy}(y^2) \right|, \qquad y>0. \]

The density of a chi-square random variable with one degree of freedom is

\[ f_X(x) = \frac{1}{\sqrt{2\pi x}} e^{-x/2}, \qquad x>0. \]

Substituting \(x=y^2\),

\[ \begin{aligned} f_Y(y) &= \frac{1}{\sqrt{2\pi y^2}} e^{-y^2/2}(2y) \\ &= \sqrt{\frac{2}{\pi}} e^{-y^2/2}, \qquad y>0. \end{aligned} \]

Thus,

\[ f_Y(y) = \sqrt{\frac{2}{\pi}} e^{-y^2/2}, \qquad y>0, \]

which is the density of a half-normal distribution^

let

\[ S= \begin{cases} 1, & \text{with probability }0.5,\\ -1, & \text{with probability }0.5, \end{cases} \]

where \(S\) is independent of \(Y\). Define

\[ Z=SY. \]

Then \(Z=\sqrt{X}\) with probability \(0.5\) and \(Z=-\sqrt{X}\) with probability \(0.5\).

For any \(z\in\mathbb{R}\),

\[ \begin{aligned} f_Z(z) &= \frac{1}{2} \sqrt{\frac{2}{\pi}} e^{-z^2/2} \\ &= \frac{1}{\sqrt{2\pi}} e^{-z^2/2}, \qquad -\infty<z<\infty. \end{aligned} \]

This is exactly the density of the standard normal distribution. Therefore,

\[ \boxed{Z\sim N(0,1).} \]

which makes sense as we are essentially flipping a coin to whether we sample from the pos/neg side of the normal dist.

2.2 (b)

set.seed(123)

n <- 10000

# Generate Chi-square(1) random variables
X <- rchisq(n, df = 1)

# Generate random signs (-1 or 1)
S <- sample(c(-1, 1),
            size = n,
            replace = TRUE,
            prob = c(0.5, 0.5)) # equal prob : .5, .5 (flip coin)

# Variable transformation
Z <- S * sqrt(X)

# Histogram
hist(Z,
     probability = TRUE,
     breaks = 40,
     main = "Standard Normal via Transformation",
     xlab = "Z",
     col = "lightblue",
     border = "white")

# Overlay the true standard normal density
curve(dnorm(x),
      from = -4,
      to = 4,
      add = TRUE,
      col = "red",
      lwd = 2)

3 Q3

3.1 (a)

Suppose

\[ X \sim \text{Gamma}(\alpha,\beta), \]

with probability density function :

\[ f_X(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}, \qquad x>0. \]

Let

\[ Y=\frac{1}{X} \]

Since

\[ x=\frac{1}{y}, \]

we have

\[ \frac{dx}{dy} = -\frac{1}{y^2}, \]

\[ \left|\frac{dx}{dy}\right| = \frac{1}{y^2}. \]

Using the change-of-variables formula :

\[ \begin{aligned} f_Y(y) &= f_X\left(\frac{1}{y}\right) \left|\frac{dx}{dy}\right| \\ &= \frac{\beta^\alpha}{\Gamma(\alpha)} \left(\frac{1}{y}\right)^{\alpha-1} e^{-\beta/y} \cdot \frac{1}{y^2}. \end{aligned} \]

Combining the powers of \(y\),

\[ \left(\frac{1}{y}\right)^{\alpha-1} \frac{1}{y^2} = \frac{1}{y^{\alpha+1}} \]

Therefore,

\[ \boxed{ f_Y(y) = \frac{\beta^\alpha}{\Gamma(\alpha)} y^{-(\alpha+1)} e^{-\beta/y}, \qquad y>0. } \]

Hence, the probability density function of \(Y=1/X\) is

\[ f_Y(y) = \frac{\beta^\alpha}{\Gamma(\alpha)} y^{-(\alpha+1)} e^{-\beta/y}, \qquad y>0. \]

3.2 (b)

set.seed(123)

n <- 10000
alpha <- 2
beta <- 3

# Target density of Y = 1/X
f <- function(y) {
  (beta^alpha / gamma(alpha)) *
    y^(-alpha - 1) *
    exp(-beta / y)
}

# Proposal density
g <- function(y) {
  1 / (1 + y)^2
}

# Maximum of f(y) / g(y)
M <- 9 * exp(-sqrt(3)) *
  (1 / sqrt(3)) *
  (1 + 1 / sqrt(3))^2

# Store accepted values of Y
y_samples <- rep(NA, n)

i <- 1

while (i <= n) {

  # Generate from g(y) = 1 / (1 + y)^2
  U1 <- runif(1)
  y_prop <- U1 / (1 - U1)

  # Accept or reject
  U2 <- runif(1)

  if (U2 <= f(y_prop) / (M * g(y_prop))) {
    y_samples[i] <- y_prop
    i <- i + 1
  }
}

# Convert Y back to X
x_samples <- 1 / y_samples

# Histogram of generated X values
hist(x_samples,
     probability = TRUE,
     breaks = 40,
     main = "Rejection Samples from Gamma(2, 3)",
     xlab = "X")

# Overlay the true Gamma density
curve(dgamma(x, shape = alpha, rate = beta),
      add = TRUE,
      col = "red",
      lwd = 2)

3.3 (b)

# Generate Gamma samples directly using rgamma()
x_rgamma <- rgamma(n,
                  shape = alpha,
                  rate = beta)

# Compare summary statistics
comparison <- data.frame(
  Method = c("Rejection method", "rgamma()"),
  Mean = c(mean(x_samples), mean(x_rgamma)),
  Variance = c(var(x_samples), var(x_rgamma)),
  Median = c(median(x_samples), median(x_rgamma))
)

comparison

##             Method      Mean  Variance    Median
## 1 Rejection method 0.6679702 0.2238089 0.5644670
## 2         rgamma() 0.6683720 0.2222580 0.5555361

# Theoretical values
theoretical_mean <- alpha / beta
theoretical_variance <- alpha / beta^2

theoretical_mean

## [1] 0.6666667

theoretical_variance

## [1] 0.2222222

as we can see our rej. method and rgamma() both are rather near and slighly different from our theoretical expectation

hist(x_samples,
     probability = TRUE,
     breaks = 40,
     xlim = range(c(x_samples, x_rgamma)),
     main = "Rejection Method vs. rgamma()",
     xlab = "X",
     col = rgb(0, 0, 1, 0.4),
     border = "white")

hist(x_rgamma,
     probability = TRUE,
     breaks = 40,
     add = TRUE,
     col = rgb(1, 0, 0, 0.4),
     border = "white")

curve(dgamma(x, shape = alpha, rate = beta),
      add = TRUE,
      col = "black",
      lwd = 2)

legend("topright",
       legend = c("Rejection method",
                  "rgamma()",
                  "True density"),
       fill = c(rgb(0, 0, 1, 0.4),
                rgb(1, 0, 0, 0.4),
                NA),
       border = c(NA, NA, NA),
       lty = c(NA, NA, 1),
       col = c(NA, NA, "black"),
       lwd = c(NA, NA, 2))