Assignment 2 Solution

Instructions: Solutions to problems 1 and 2 are to be submitted on Quercus (PDF les 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 Z 1




by taking a bounded function h with h(x) 0 for all x and h(x) dx < 1. We then de ne the region




 

q




Ch = (u; v) : 0 u h(v=u)




and generate (U; V ) uniformly distributed on Ch. We then de ne the random variable X = V=U.




(a) The joint density of (U; V ) is




f(u; v) =
1
for (u; v) 2 Ch
jChj



where jChj is the area of Ch. Show that the joint density of (U; X) is








u


for 0 u q


g(u; x) =




h(x)


h




jC
j





and that the density of X is h(x) for some 0.




The implementation of this method is somewhat complicated by the fact that it is typically di cult to sample (U; V ) from a uniform distribution on Ch. However, it is usually



possible to nd a rectangle of the form Dh = f(u; v) : 0 u u+; v v v+g 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 de ne u+, v and
v+ as follows:




























+ =
x
q


x
q




+
x
q




u


h(x) v
h(x):


max


h(x) v = min x




= max x





(Hint: It su ces to show that if (u; v) 2 Ch then (u; v) 2 Dh where Dh is de ned using u+,

, 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 =
q


= 0:8577639, and v+ = q




2=e
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 f"ig is a sequence of random variables with mean 0 and nite variance representing noise. We will assume that 1; ; n are smooth in the sense that i = g(i) for some continuous and di erentiable function g. The least squares estimates of 1; ; n are trivial




bi = yi for all i | but we can modify least squares in a number of ways to accommodate the \smoothness" assumption on f ig. In this problem, we will consider estimating f ig by
minimizing

n 1
(yi i)2
+
( i+1 2 i + i 1)2
X


Xi
i=1


=2



where 0 is a tuning parameter that controls the smoothness of the estimates b1; ; bn. (This method is known as Whittaker graduation in actuarial science and the Hodrick-Prescott lter in economics.)




Show that if fyig are exactly linear, i.e. yi = a i + b for all i then bi = yi for all i.



In principal, we could compute b1; ; bn using ordinary least squares estimation. Show
that b = (b1; ; bn)T minimizes

ky X k2




where y is a vector of length 2n 2 and X is an (2n 2)


n matrix. What are y and X?





When n is large, computing b1; ; bn 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 modi cation of the Gauss-Seidel algorithm, which at



each stage solves a p( n) variable least squares problem. The basic algorithm is as follows:




Initialize b.



Randomly sample a subset w of size p from the integers 1; ; n. De ne 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; de ne w and w analogously so that X =



Xw w + Xw w.




De ne bw to minimize


Xw
w
Xw w


y






b




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.




On Quercus, there is a function HP in a le HP.txt, which implements the randomized block Gauss-Seidel algorithm outlined in part (c). This function can be used as follows



r <- HP(x,lambda=2000,p=20,niter=100)




where lambda is the value of , p is the value of p, and niter speci es the number of iterations of the algorithm. The output of the function (contained here in r) consists of two components: the estimates of f g 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.




(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 nd 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 oating point operations needed increases with p like p2.


















































































3
Supplemental problems (not to hand in):







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 nite 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 nite set is f0; 1; ; n 1g for some integer n 2. For U Unif(0; 1), de ne




Y = bnUc and V = nU Y




where bxc is the integer part of x. Show that Y has a uniform distribution on the set f0; 1; ; n 1g, V has a uniform distribution on [0; 1], and Y and V are independent.

(b) What happens to the precision of V de ned in part (a) as n increases? (For example, if




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?






One issue with rejection sampling is a lack of e ciency 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:




Generate two independent random variables U Unif(0; 1) and X with density g.



If U f1(X)=g(X) then return X.



Otherwise (that is, if U f1(X)=g(X)) generate X from the density



f2 (x) =
f2(x)
:




Z 1 f2(t) dt



Note that we must be able to easily sample from g and f2 in order for the A-C method to be e cient; in some cases, they can both be taken to be uniform distributions.




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


( 1 x
1)












(1 + x2)




2








x
2
using the A-C method with f (x) = k, a constant, for 1




1 (so that f (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) De ning 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?







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 constant b. The density in question is
exp( x2=2)


f(x) = p2 (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.)




De ne 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 de nes X =



b + Y =b if

Y 2

2 ln(U) b2 :







(Hint: Note that b + Y =b has density g.)




(b) Show the probability of acceptance is given by

p




2 b(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 de ned 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 max f(x)




0 x b 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 min
f(x)












g (x)












x b 0








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






















kxk
1 exp(
x)
for x 0.
gk(x) =








(k 1)!
Show that






k
1 j exp( )








P (Tk 1) = Z11 gk(x) dx = j=0
j!








X












(Hint: Use integration by parts.)




How might you use the result of part (a) to generate random variables from a Poisson distribution with mean ?




























































6