Question

R Code

2. Suppose X is a random variable with density f(x; ψ) = 2ψ-2 xI(0 ≤ x ≤ ψ), where ψ > 0 is a parameter.

(a) Draw a graph of the density when ψ = 3 in R.

(b) Write a function called genTri that generates n independent realizations of X. The function should take the arguments n and psi and return a vector of n realizations. Make a histogram of n = 104 realizations and compare the histogram to the density plot in the previous question.

(c) Derive an expression for E(X) – it will be a function of ψ. Confirm your expression when ψ = 3 using a Monte Carlo simulation.

(d) A 95% confidence interval for E(X) is given by Xbar ± t0.975,n−1 * S/√ n, where Xbar is the sample mean of n independent copies of X and t0.975,n−1 denotes the 97.5th percentile of the t-distribution with n − 1 degrees of freedom. For n = 5, 10, 15, 100 and ψ = 3, estimate the probability P(Xbar − t0.975,n−1 * S/√ n ≤ ψ ≤ Xbar + t0.975,n−1 * S/√ n) using a Monte Carlo simulation with 104 replications. Note: you can use the qt function in R to generate the respective t percentiles.

Answer 1

The pdf of the random variable is

a) For , the graph of the pdf is plotted below using R.

curve(2*x/9, ylab = "f(x)", xlim = c(0,3), col="darkred", lwd=2)

b) Using Invesre transform method, we can generate random variables from the above distribution.

If is standard uniform RV, then is RV from .

The R code (with function genTri) to generate 104 realizations of the RV and plot the histogram is given below.

n_random = 104
psi <- 3
genTri <- function(n, psi)
{
U <- runif(n)
X <- sqrt(U)*psi
return(X)
}
X <- genTri(n_random, psi)
hist(X, prob = TRUE,freq = FALSE, col = "skyblue", xlab = "x", main = "Histogram")

The histogram has increasing steps.

c) The expected value of the RV is

When , theoretical mean is and the simulated mean is .

Codee below.

n_random = 104
psi <- 3
genTri <- function(n, psi)
{
U <- runif(n)
X <- sqrt(U)*psi
return(X)
}
X <- genTri(n_random, psi)
mean(X)

d) The R code for finding the probability

where is given below.

nsims = 104
n <- 100
psi <- 3
n_total <- 0
genTri <- function(n, psi)
{
U <- runif(n)
X <- sqrt(U)*psi
return(X)
}
for(i in 1:nsims)
{
X <- genTri(n, psi)
xbar <- mean(X)
s <- sd(X)
t <- qt(0.975,n-1)
if (xbar - t*s/sqrt(n) < 2*psi/3 & xbar + t*s/sqrt(n) > 2*psi/3 )
{
n_total <- n_total + 1
}
  
}

n_total/nsims

The probability for

is

is

is

is