Description
1. In this question, you will create a rejection sampling algorithm to generate a realization of a random variable with density f defined by f(x) = 3 4 (1 − x 2 ) if x ∈ (−1, 1) 0 otherwise . (a) Create this rejection sampling algorithm using Unif(−1, 1) as the trial/proposal distribution. Simplify and state the steps of this algorithm. (b) Write an R function called rquad that generates a realization of n independent copies of the random variable X with density f. This function must use the algorithm created in part 1a. This function has one argument n, which is the sample size. This function returns a list with two elements: • x.list, a vector of n entries, where the ith entry is the realization of the ith independent copy of X. • k.list, a vector of n entries, where the ith entry is the number of iterations of the rejection sampling algorithm required to produce the realization of the ith independent copy of X. (c) Test rquad by generating a realization of X1, . . . , X1000 iid with density f. Create a a histogram of these measurements. How many iterations, on average, were required to produce each realization? 2. (a) Write an R function called myrtnorm that generates a realization of the sequence of independent random variables T1, . . . , Tn, where each Ti has the truncated Normal distribution: Ti ∼ (X|a < X < b) where X ∼ N(µ, σ2 ) and a, b are non-random. A rejection sampling algorithm should be used in combination with the Box–Muller method. Only the standard uniform random generator runif is allowed. This function should have five arguments: • n, the random sample size • mu, the value of µ • sigma, the value of σ • a, the value of a (the left endpoint of the interval or -Inf) • b, the value of b (the right endpoint of the interval or Inf). This function should return a vector of n entries with the generated realization of T1, . . . , Tn. (b) Derive an expression for the expected number of iterations of the rejection sampling algorithm required to generate a realization of T1 in terms of µ, σ, a and b. 1 (c) Fix n = 500 and pick values of µ, σ, a and b so that the expected number of iterations of the rejection sampling algorithm required to generate a realization of T1 is less than 10. Using these values, use myrtnorm to generate a realization of T1, . . . , T500 and produce a histogram. 3. (a) Write an R function called myrexp that generates a realization of n independent random variables X1, . . . , Xn, where Xi has the exponential distribution with mean µ > 0 for i = 1, . . . , n. In other words, generate a realization of a random sample of size n from Exp(µ). Only calls to R’s standard uniform generator runif are permitted, e.g. calling rexp is not allowed. This function should have two arguments: • n, the random sample size • mu, the user-specified mean of the exponential distribution This function should return a vector of n entries with the generated realization of X1, . . . , Xn. (b) Test myrexp by generating a realization of random sample of size 1000 from the Exponential distribution with some mean µ that you pick. Create a QQ–plot to compare the data percentiles (of the realization of the random sample) to the percentiles of Exp(µ). Calling the function qexp is not allowed here. (c) Let Y1, . . . , Yn be independent copies of Y ∼ Exp(µ) and define Y¯ = n −1 Pn i=1 Yi . Write an R function called run.exp.sim that generates a realization of reps independent copies of Y¯ with sample size n. This function should have three arguments: • n the sample size • mu the mean of the exponential distribution • reps the number of realizations of Y¯ The function should return a vector of reps realizations of Y¯ and display a Normal QQ–plot of entries in this vector. (d) A civil engineer measured the times between vehicle arrivals at a rural bridge on a Sunday afternoon. For a simple model, she assumes that her measured inter-arrival times (in minutes) x1, . . . , x30 are a realization of a random sample from the exponential distribution with unknown mean µ. She computes the observed sample mean ¯x = (1/30)P30 i=1 xi to estimate µ. Is this sample size of 30 large enough for ¯x to be a realization of random variable with a distribution well approximated by the Normal distribution? To respond, pretend that µ = 3.4 minutes and perform a simulation study using the function run.exp.sim. Comment on the result. 4. (a) Write an R function called mymvrnorm that generates a realization of the sequence of independent random vectors Y1, . . . , Yn, where Yi ∼ Np(µ, Σ) for i = 1, . . . , n. Only calls to R’s standard uniform generator runif are permitted. This function should have three arguments: • n, the random sample size • mu, the mean vector with p entries 2 • Sigma, this is the covariance matrix Σ ∈ S p 0 . This function should return a matrix with n rows and p columns, where the ith row has the realization of Yi . The eigen function should be called in your definition of mymvrnorm. (b) Suppose that we plan to measure the heights of n individuals. i. Suppose that the yet-to-be measured heights X1, . . . , Xn are a random sample from N(µ, σ2 ). Compute the mean and variance of X¯ = n −1 Pn i=1 Xi (the standard estimator of µ). Express these in terms of µ, σ and n. ii. Suppose that the yet-to-be measured heights will be of individuals in the same family. Let (H1, . . . , Hn) 0 be these yet-to-be measured heights and suppose they have n-variate normal distribution with mean vector (µ, . . . , µ) 0 ∈ R n and covariance matrix Σ ∈ S n 0 with (i, j)th entry Σij = σ 2 · 0.7 |i−j| for (i, j) ∈ {1, . . . , n} × {1, . . . , n}. Compute the mean and variance of H¯ = n −1 Pn i=1 Hi (the standard estimator of µ). Express these in terms of µ, σ and n. iii. Is H¯ better or worse than X¯ as an estimator of µ. Explain. (c) Design and perform a simulation study that compares the formulas for the mean and variance of X¯ and H¯ derived in part 4b to their corresponding simulated estimates. You should use mymvrnorm. Set n = 10, µ = 68, and σ = 3. 3
