Assignment #2 STA410H1F/2102H1F due Friday October 18, 2019 Instructions: Solutions to problems 1 and 2 are to be submitted on Quercus (PDF files only). You are strongly encouraged to do problems 3–6 but these are not to be submitted for grading. 1. An interesting variation of rejection sampling is the ratio of uniforms method. We start by taking a bounded function h with h(x) ≥ 0 for all x and ∫ ∞ −∞ h(x) dx < ∞. We then define the region Ch = { (u, v) : 0 ≤ u ≤ √ h(v/u) } and generate (U, V ) uniformly distributed on Ch. We then define the random variable X = V/U . (a) The joint density of (U, V ) is f(u, v) = 1 |Ch| for (u, v) ∈ Ch where |Ch| is the area of Ch. Show that the joint density of (U,X) is g(u, x) = u |Ch| for 0 ≤ u ≤ √ h(x) and that the density of X is γ h(x) for some γ > 0. (b) The implementation of this method is somewhat complicated by the fact that it is typically difficult to sample (U, V ) from a uniform distribution on Ch. However, it is usually possible to find a rectangle of the form Dh = {(u, v) : 0 ≤ u ≤ u+, v− ≤ v ≤ v+} such that Ch is contained within Dh. Thus to draw (U, V ) from a uniform distribution on Ch, we can use rejection sampling where we draw proposals (U∗, V ∗) from a uniform distribution on the rectangle Dh; note that the proposals U∗ and V ∗ are independent random variables with Unif(0, u+) and Unif(v−, v+) distributions, respectively. Show that we can define u+, v− and v+ as follows: u+ = max x √ h(x) v− = min x x √ h(x) v+ = max x x √ h(x). (Hint: It suffices to show that if (u, v) ∈ Ch then (u, v) ∈ Dh where Dh is defined using u+, v−, and v+ above.) (c) Implement (in R) the method above for the standard normal distribution taking h(x) = exp(−x2/2). In this case, u+ = 1, v− = − √ 2/e = −0.8577639, and v+ = √ 2/e = 0.8577639. What is the probability that proposals are accepted? 1 2. Suppose we observe y1, · · · , yn where yi = θi + εi (i = 1, · · · , n) where {εi} is a sequence of random variables with mean 0 and finite variance representing noise. We will assume that θ1, · · · , θn are smooth in the sense that θi = g(i) for some continuous and differentiable function g. The least squares estimates of θ1, · · · , θn are trivial — θ̂i = yi for all i — but we can modify least squares in a number of ways to accommodate the “smoothness” assumption on {θi}. In this problem, we will consider estimating {θi} by minimizing n∑ i=1 (yi − θi)2 + λ n−1∑ i=2 (θi+1 − 2θi + θi−1)2 where λ > 0 is a tuning parameter that controls the smoothness of the estimates θ̂1, · · · , θ̂n. (This method is known as Whittaker graduation in actuarial science and the Hodrick- Prescott filter in economics.) (a) Show that if {yi} are exactly linear, i.e. yi = a× i+ b for all i then θ̂i = yi for all i. (b) In principal, we could compute θ̂1, · · · , θ̂n using ordinary least squares estimation. Show that θ̂ = (θ̂1, · · · , θ̂n)T minimizes ‖y∗ −Xθ‖2 where y∗ is a vector of length 2n−2 and X is an (2n−2)×n matrix. What are y∗ and X? (c) When n is large, computing θ̂1, · · · , θ̂n directly, for example, using the OLS formulation in part (b) is computationally expensive when n is large. Alternatively, we could use the Gauss-Seidel algorithm but it converges slowly, particularly for larger values of λ. One possible alternatively is a randomized modification of the Gauss-Seidel algorithm, which at each stage solves a p( n) variable least squares problem. The basic algorithm is as follows: 0. Initialize θ̂. 1. Randomly sample a subset w of size p from the integers 1, · · · , n. Define Xw to be the submatrix of X with column indices w and Xw¯ to be the submatrix of X with column indices in the complement of w; define θw and θw¯ analogously so that Xθ = Xwθw +Xw¯θw¯. 2. Define θ̂w to minimize ∥∥∥y∗ −Xw¯θ̂w¯ −Xwθw∥∥∥ with respect to θw. 2 3. Repeat steps 1 and 2 until convergence. Show that the objective function is non-increasing from one iteration to the next. (d) On Quercus, there is a function HP in a file HP.txt, which implements the randomized block Gauss-Seidel algorithm outlined in part (c). This function can be used as follows > r where lambda is the value of λ, p is the value of p, and niter specifies the number of iterations of the algorithm. The output of the function (contained here in r) consists of two components: the estimates of {θ} in r$xhat and the values of the objective function for each iteration in r$ss. Use this function (with λ = 2000) on data on monthly yields on short-term British securities over a 21 year period (252 months). Try out various values of p between 5 and 50 (using 1000 iterations). Describe how the objective function decreases with each iteration as p varies between 5 and 50. (e) (Optional) Methods such as the randomized block Gauss-Seidel algorithm are essential in problems where the number of unknown parameters is extremely large. The goal is not to minimize the objective function but merely to find a solution that’s close to the minimum. In the context of the randomized block Gauss-Seidel algorithm, what factors should you consider in choosing p, the numbers of parameters being optimized at each step? Keep in mind that for least squares, the number of floating point operations needed increases with p like p2. 3 Supplemental problems (not to hand in): 3. To generate random variables from some distributions, we need to generate two indepen- dent two independent random variables Y and V where Y has a uniform distribution over some finite set and V has a uniform distribution on [0, 1]. It turns out that Y and V can be generated from a single Unif(0, 1) random variable U . (a) Suppose for simplicity that the finite set is {0, 1, · · · , n− 1} for some integer n ≥ 2. For U ∼ Unif(0, 1), define Y = bnUc and V = nU − Y where bxc is the integer part of x. Show that Y has a uniform distribution on the set {0, 1, · · · , n− 1}, V has a uniform distribution on [0, 1], and Y and V are independent. (b) What happens to the precision of V defined in part (a) as n increases? (For example, if U has 16 decimal digits and n = 106, how many decimal digits will V have?) Is the method in part (a) particularly feasible if n is very large? 4. One issue with rejection sampling is a lack of efficiency due to the rejection of random variables generated from the proposal density. An alternative is the acceptance-complement (A-C) method described below. Suppose we want to generate a continuous random variable from a density f(x) and that f(x) = f1(x) + f2(x) (where both f1 and f2 are non-negative) where f1(x) ≤ g(x) for some density function g. Then the A-C method works as follows: 1. Generate two independent random variables U ∼ Unif(0, 1) and X with density g. 2. If U ≤ f1(X)/g(X) then return X. 3. Otherwise (that is, if U > f1(X)/g(X)) generate X from the density f ∗2 (x) = f2(x)∫ ∞ −∞ f2(t) dt . Note that we must be able to easily sample from g and f ∗2 in order for the A-C method to be efficient; in some cases, they can both be taken to be uniform distributions. (a) Show that the A-C method generates a random variable with a density f . What is the probability that the X generated in step 1 (from g) is “rejected” in step 2? 4 (b) Suppose we want to sample from the truncated Cauchy density f(x) = 2 pi(1 + x2) (−1 ≤ x ≤ 1) using the A-C method with f2(x) = k, a constant, for −1 ≤ x ≤ 1 (so that f ∗2 (x) = 1/2 is a uniform density on [−1, 1]) with f1(x) = f(x)− f2(x) = f(x)− k (−1 ≤ x ≤ 1). If g(x) is also a uniform density on [−1, 1] for what range of values of k can the A-C method be applied? (c) Defining f1, f2, and g as in part (b), what value of k minimizes the probability that X generated in step 1 of the A-C algorithm is rejected? 5. Suppose we want to generate a random variable X from the tail of a standard normal distribution, that is, a normal distribution conditioned to be greater than some constan
t b. The density in question is f(x) = exp(−x2/2)√ 2pi(1− Φ(b)) for x ≥ b with f(x) = 0 for x < b where Φ(x) is the standard normal distribution function. Consider rejection sampling using the shifted exponential proposal density g(x) = b exp(−b(x− b)) for x ≥ b. (This proposal density is used by the Monty Python algorithm to sample from the tail of the normal distribution.) (a) Define Y be an exponential random variable with mean 1 and U be a uniform random variable on [0, 1] independent of Y . Show that the rejection sampling scheme defines X = b+ Y/b if −2 ln(U) ≥ Y 2 b2 . (Hint: Note that b+ Y/b has density g.) (b) Show the probability of acceptance is given by √ 2pib(1− Φ(b)) exp(−b2/2) . What happens to this probability for large values of b? (Hint: You need to evaluate M = max f(x)/g(x).) 5 (c) Suppose we replace the proposal density g defined above by gλ(x) = λ exp(−λ(x− b)) for x ≥ b. (Note that gλ is also a shifted exponential density.) What value of λ maximizes the proba- bility of acceptance? (Hint: Note that you are trying to solve the problem min λ>0 max x≥b f(x) gλ(x) for λ. Because the density gλ(x) has heavier tails, the minimax problem above will have the same solution as the maximin problem max x≥b min λ>0 f(x) gλ(x) which may be easier to solve.) 6. (a) Suppose that E1, E2, · · · are independent Exponential random variables with density f(x) = λ exp(−λx) for x ≥ 0. Then the distribution of Tk = E1 + · · · + Ek is a Gamma distribution whose density is gk(x) = λkxk−1 exp(−λx) (k − 1)! for x ≥ 0. Show that P (Tk ≥ 1) = ∫ ∞ 1 gk(x) dx = k−1∑ j=0 λj exp(−λ) j! (Hint: Use integration by parts.) (b) How might you use the result of part (a) to generate random variables from a Poisson distribution with mean λ? 6
