Sampling From Discrete Distributions

1 Sampling from discrete distributions A discrete random variable X is a random variable that has a probability mass function p(x) = P (X = x) for any x ∈ S, where S = {x1 , x2 , ..., xk } denotes the sample space, and k is the (possibly infinite) number of possible outcomes for the discrete variable X, and suppose S is ordered from smaller to larger values. Then the CDF, F for X is X F (xj ) = p(xi ) i≤j Discrete random variables can be generated by slicing up the interval (0, 1) into subintervals which define a partition of (0, 1): (0, F (x1 )), (F (x1 ), F (x2 )), (F (x2 ), F (x3 )), ..., (F (xk−1 ), 1), generating U = Uniform(0, 1) random variables, and seeing which subinterval U falls into. Let Ij = I (U ∈ (F (xj ), F (xj+1 ))). Then,   P (Ij = 1) = P F (xj−1 ) ≤ U ≤ F (xj ) = F (xj ) − F (xj−1 ) = p(xj ) where F (x0 ) is defined to be 0. So, the probability that Ij = 1 is same as the probability that X = xj , and this can be used to generate from the distribution of X. As an example, suppose that X takes values in S = {1, 2, 3} with probability mass function defined by the following table: p(x) p1 p2 p3 x 1 2 3 To generate from this distribution we partition (0, 1) into the three sub-intervals (0, p1 ), (p1 , p1 + p2 ), and (p1 + p2 , p1 + p2 + p3 ), generate a Uniform(0, 1), and check which interval the variable falls into. The following R code does this, and checks the results for p1 = .4, p2 = .25, and p3 = .35: # n is the sample size # p is a 3-length vector containing the corresponding probabilities rX <- function(n, p) { # generate the underlying uniform(0,1) variables U <- runif(n) # storage for the random variables, X X <- rep(0,n) # which U are in the first interval w1 <- which(U <= p[1]) X[w1] <- 1 # which U are in the second interval w2 <- which( (U > p[1]) & (U < sum(p[1:2])) ) X[w2] <- 2 # which U are in the third interval w3 <- which( U > sum(p[1:2]) ) X[w3] <- 3 return(X) } X = rX(10000, c(.4, .25, .35)) mean(X == 1) [1] 0.407 mean(X == 2) [1] 0.2504 mean(X == 3) [1] 0.3426 The empirical probabilities appear to agree with the true values. Any discrete random variable with a finite sample space can be generated analogously, although the use of a for loop will be necessary when the number of intervals to check is large. The process will be similar when the variable has an infinite sample space– one example of this is the Poisson distribution. The probability mass function for the poisson with parameter λ has the form e−λ λx p(x) = x! whose sample space is all non-negative integers. The following R program generates from this distribution and compared the empirical mass function with the true mass function for λ = 4: # function to calculate poission CDF at an integer x # when lambda = L. This does the same thing as ppois() pois.cdf <- function(x, L) { # grid of values at which to calculate the mass function v <- seq(0, x, by=1) # return CDF value return( sum( exp(-L)*(L^v)/factorial(v) ) ) } # n: sample size, L: lambda r.pois <- function(n, L) { U <- runif(n) X <- rep(0,n) # loop through each uniform for(i in 1:n) { # first check if you are in the first interval if(U[i] < pois.cdf(0,L)) { X[i] <- 0 } else { # while loop to determine which subinterval,I, you are in # terminated when B = TRUE B = FALSE I = 0 while(B == FALSE) { # the interval to check int <- c( pois.cdf(I, L), pois.cdf(I+1,L) ) # see if the uniform is in that interval if( (U[i] > int[1]) & (U[i] < int[2]) ) { # if so, quit the while loop and store the value X[i] <- I+1 B = TRUE } else { # If not, continue the while loop and increase I by 1 I=I+1 } } } } return(X) } # generate 1000 Pois(4) random variables V = r.pois(1000, 4) # empirical mass function Mf <- c(0:15) for(i in 0:15) Mf[i+1] <- mean(V==i) # plot the observed mass function with the # the true mass function overlaying it b <- c(0:15) plot(b, Mf, xlab="x", ylab="p(x)", main="Empirical mass function", col=2) lines(b, Mf, col=2) # overlay with true mass function points(b, dpois(b,4), col=4) lines(b, dpois(b,4), col=4) 2 Sampling from continuous distributions Continuous random variables are (informally) those whose sample space is composed of real intervals not exclusively containing integers. In a continuous distribution the probability of taking on any particular value in the sample space in 0; probabilities can only be assigned to intervals in a continuous distribution. For this reason the logic of the previous section does not apply directly and other methods must be used. 2.1 The inversion method It is a fact that if X has CDF F , then F (X) has a Uniform(0, 1) distribution. The proof of this is a straightforward calculation: P (F (X) ≤ x) = P (F −1 (F (X)) ≤ F −1 (x)) = P (X ≤ F −1 (x)) = F (F −1 (x)) =x So the CDF of F (X) is x, which is the same as the CDF of as Uniform(0, 1). Here F −1 denotes the inverse of the CDF (also called the quantile function) and is defined as the function which satisfies: F (x) = y ⇐⇒ x = F −1 (y) In other words, F (F −1 (y)) = y and F −1 (F (x)) = x. The fact above implies that if X has CDF F , then F −1 (U ) will have CDF F . So, if you are able to calculate F −1 , and can generate uniforms, then you can generate a sample from F . As an example, suppose you want to generate from a distribution with CDF x 1+x for x ∈ (0, ∞). This distribution has mean and variance equal to ∞. To calculate F −1 (y), you specify a value for y and solve for x: F (x) = x 1+x (1 + x)y = x y = x − xy y = x(1 − y) y =x y−1 y= so F −1 (y) = y . y−1 Therefore U 1−U # n is the sample size r.harmonic <- function(n) { # generate uniforms U <- runif(n) # return F^-1(U) return( U/(1-U) ) } will have CDF F : X <- r.harmonic(1000) # empirical CDF v <- seq(0, 20, length=1000) emp.cdf <- rep(0, 1000) for(i in 1:1000) emp.cdf[i] <- mean(X <= v[i]) # true CDF true.cdf <- v/(1+v) plot(v, emp.cdf, xlab="X", ylab="F(X)", main="Empirical vs True CDF", col=2,type="l") lines(v,true.cdf,col=4) Due to the long tails, this distribution is a good candidate for a trial distribution in rejection sampling, which we will mention later. As a second example suppose X has CDF  F (x) = 1 1 + e−x θ where θ > 0 is a parameter. This distribution is known as the skew logistic distribution, which is symmetric when θ = 1, and skewed otherwise. This time the inversion is a bit more involved: θ 1 y= 1 + e−x 1/y = (1 + e−x )θ  (1/y)1/θ = 1 + e−x
− log (1/y)1/θ − 1 = x
So F −1 (y) = − log (1/y)1/θ − 1 . # generate n samples from the skew logistic # distribution with parameter Theta r.skewlog <- function(n, Theta) { U <- runif(n) return( -log( (1/U)^(1/Theta) - 1) ) } X <- r.skewlog(1000, 4) # empirical CDF v <- seq(-10, 10, length=1000) emp.cdf <- rep(0, 1000) for(i in 1:1000) emp.cdf[i] <- mean(X <= v[i]) # true CDF true.cdf <- (1 + exp(-v))^(-4) plot(v, emp.cdf, xlab="X", ylab="F(X)", main="Empirical vs True CDF", col=2,type="l") lines(v,true.cdf,col=4) In many cases you can not symbolically invert the CDF (the normal distribution is an example of this). In principle, you can still use this method in such situations, but you will have to numerically calculate the quantile function. When we get to the section of the course dealing with line searching (optimization), we will do an example like this. 2.2 Transformation method In some situations you can know mathematically that a particular function of a random variable has a certain distribution. On homework 1, problem 2 you were given an example of this– the transformation introduced there is called the Box-Muller transform, and is actually the preferred way to generate normally distributed variables. When such relationships are know, it gives a simple way of generating from a distribution. In this example, suppose we wish to generate from the exponential(θ) distribution, and only have access to a computer which generates numbers from the skew logistic distribution, and do not know the inversion method. It turns out that if X ∼ SkewLogistic(θ), then log(1 + e−X ) is exponential(θ): P (log(1 + e−X ) ≤ y) = P (1 + e−X ≤ ey ) = P (e−X ≤ ey − 1) = P (−X ≤ log(ey − 1)) = P (X ≥ − log(ey − 1)) = 1 − P (X ≤ − log(ey − 1))  θ 1 =1− 1 + elog(ey −1)  θ 1 =1− 1 + ey − 1 = 1 − e−yθ which is the CDF of an exponential(θ). We can check by simulation that this transformation is correct: X <- r.skewlog(1000, 3) Y <- log(1 + exp(-X)) Several other useful transformation identities are: P • if X1 , ..., Xk ∼ N (0, 1), then ki=1 Xi2 ∼ χ2k . P • if X1 , ..., Xk ∼ exponential(λ), then ki=1 Xi2 ∼ Gamma(k, λ) • if U ∼ Uniform(0, 1), then 1 − X 1/k ∼ Beta(1, n) • if X1 ∼ Gamma(α, λ) and X2 ∼ Gamma(β, λ), then X1 X2 +X1 ∼ Beta(α, β) • if X ∼ t(df ), then X 2 ∼ F (1, df ) As an exercise for later check by simulation that these identities hold– some useful R functions will be rnorm(), rnorm(), rchisq(), rbeta(), rgamma(), rf(), rt(). 2.3 Rejection Sampling Rejection sampling is a general algorithm to generate samples from a distribution with density p1 (called the target density) based only on the ability to generate from another distribution with density p2 (usually referred to as the trial density) such that sup x p1 (x) ≤M <∞ p2 (x) This is essentially saying that the tails of the trial density, p2 , are not infinitely fatter than those of the target density, p1 . The basic algorithm is: 1. Generate U ∼ Uniform(0, 1) 2. Generate X ∼ p2 3. If U ≤ Mp·p1 (X) then accept X as a realization from p1 , otherwise throw out X and 2 (X) try again Showing the distribution of accepted draws is the same as the target distribution is a simple conditional probability calculation. First we will show that the rejection rate is 1/M : P (accept) = X P (Accept|X = x)P (Drawing X = x) x = X  P x = X x =  p1 (x) |X = x p2 (x) U≤ M · p2 (x) p1 (x) · p2 (x) M · p2 (x) X p1 (x) x M = 1/M where P (accept|X) = Mp1p(X) , follows from the basics about the uniform CDF. Now to 2 (X) show that the distribution of accepted draws is the same as the target distribution: P (X, accept) P (accept) P (accept|X)P (Drawing X) = P (accept) P (accept|X)p2 (X) = 1/M p1 (X) = · M p2 (X) M p2 (X) = p1 (X) P (X|accept) = The second lines follows from Bayes’ Rule. The third line follows from the fact that P (Drawing X) = p2 (X), and P (accept) = 1/M . One drawback of this method is that you will end up generating more uniforms than necessary, since some of them will get rejected. Theoretically, as long as you know that p1 indeed has longer tails then p2 , you can choose M to be ridiculously large and this will still yield a valid algorithm. Since the rejection rate is equal to 1/M , you will have to, on average, generate M · n draws from the trial distribution and from the uniform distribution just to get n draws from the target distribution. So choosing M to be too large will yield an inefficient algorithm, and it pays to find a sharp bound for M . The Cauchy distribution is one that will have larger tails than almost any other distribution you will ever deal with, so usually makes a good trial density. Also, for most nonnegative random variables, the first distribution we demonstrated the inversion method can be used as a trial density, which we will use here. As an example we will generate from the exponential(λ) distibution using rejection sampling. Recall the CDF F (x) = x/(1 + x) from a previous example; the density corresponding to this CDF is x ∂ ∂F (x) = 1+x ∂x ∂x (1 + x) − x = (1 + x)2  2 1 = 1+x We will take this to be p2 . In this case p1 is the exponential density, p1 (x) = λe−λx . To get a sharp bound for M , we now have to determine the maximum of G(x) = p1 (x) = λe−λx (1 + x)2 p2 (x) We calculate the derivative of this function:  0 −λx −λx 2  − λe (1 + x) + 2(1 + x)e   −λx 2 = λe − λ(1 + x) + 2(1 + x)   −λx 2 = λe − λx + 2(1 − λ)x + (2 − λ) G (x) = λ You can check that when G0 (x) = 0 has a solution, it is indeed a maximum. Since λe−λx is strictly positive, the solution to G0 (x) = 0 must be the same as the solution to Q(x) = −λx2 + 2(1 − λ)x + (2 − λ) = 0, a quadratic equation with discriminant 4(1 − λ)2 − 4λ(2 − λ) which is always 4 independently of λ. So any solution of this quadratic will have the form −2(1 − λ) ± 2 (1 − λ) ± 1 = = (−1 + 2/λ, −1) −2λ λ since x ≥ 0, the only feasible solution is x? = −1 + 2/λ. Thus when x? = −1 + 2/λ ≥ 0, M = G(x? ) constitutes a sharp bound. When −1 + 2/λ ≤ 0 (i.e. λ ≤ 2), there is no feasible solution. However, when λ ≤ 2, Q(0) = (2 − λ) ≤ 0, and furthermore, Q0 (x) = −2λx + 2(1 − λ) = 2 − 2λ(1 + x) ≤ 2 − 2(1 + x) ≤0 (since λ ≤ 2) (since − 2(1 + x) ↓) and so Q is monotonically decreasing which, combined with Q(0) ≤ 0, implies Q is negative for any x. Therefore, G0 (x) = λe−λx Q(x) ≤ 0, so G is also monotonically decreasing, and x? = 0. Finally this implies M = G(0) when λ ≥ 2. To summarize: • If λ ≤ 2, M = G(−1 + 2/λ) • If λ ≥ 2, M = G(0). We can plot M as a function of λ: Lambda <- seq(.1, 5, length=1000) M <- rep(0, 1000) for(j in 1:1000) { G <- function(x) Lambda[j]*( (1+x)^2 )*exp(-Lambda[j]*x) if( Lambda[j] < 2 ) { M[j] <- G(-1 + 2/Lambda[j]) } else { M[j] <- G(0) } } plot(Lambda, M) Now that we have an explicit bound we can carry out rejection sampling. We will use λ = 3/2 in this example. We will also keep track of the total number of iterations required to collect N = 1000 samples, stored in the variable niter. Lambda = 1.5 # sample size N = 1000 # target density p1 <- function(x) Lambda*exp(-x*Lambda) # trial density p2 <- function(x) (1 + x)^(-2) # function to generate from trial density r.trial <- function(n) { U <- runif(n) return( U/(1-U) ) } # bound on p1(x)/p2(x) if(Lambda < 2) { M = p1(-1 + 2/Lambda)/p2(-1 + 2/Lambda) } else { M = p1(0)/p2(0) } M = 100 X <- rep(0, N) # number of iterations used niter = 0 for(j in 1:N) { # keep going until you get an accept B <- FALSE while(B == FALSE) { # add to the number of iterations niter = niter + 1 # generate a candidate V <- r.trial(1) # calculate the ratio, evaluated at V R <- p1(V)/(M*p2(V)) # generate a single uniform U <- runif(1) # if the ratio is lower the uniform, # accept this as a draw if(U < R) { X[j] <- V B <- TRUE } } } M [1] 1.617415 mean(X) [1] 0.6927637 var(X) [1] 0.449128 niter [1] 1625 Notice the sample mean and variance appear to agree with the theoretical mean and variance, which are 2/3, and 4/9, respectively. Also note the number of iterations required was 1625, which is approximately equal to n · M = 1617.415. Now consider if we were lazy and just chose M = 100, and carried out the sampling (cut and paste the same code but uncomment the part M=100). > mean(X) [1] 0.6609276 > var(X) [1] 0.4499482 > niter [1] 98168 Notice we still appear to have sampled from the right distribution, but we used nearly 90 times the number of iterations, emphasizing the need for a sharp bound on M .