1 Sampling from discrete distributions

1 1 Sampling from discrete distributionsA discrete random variableXis a random variable that has a probability mass functionp(x) =P(X=x) for anyx S, whereS={x1, x2, .., xk}denotes the sample space, andkis the (possibly infinite) number of possible outcomes for the discrete variableX, andsupposeSis ordered from smaller to larger values. Then the CDF,FforXisF(xj) = i jp(xi) discrete random variables can be generated by slicing up the interval (0,1) into subinter-vals which define apartitionof (0,1):(0, F(x1)),(F(x1), F(x2)),(F(x2), F(x3)).

2 ,(F(xk 1),1),generatingU= Uniform(0,1) random variables, and seeing which subintervalUfalls (U (F(xj), F(xj+1))). Then,P(Ij= 1) =P(F(xj 1) U F(xj))=F(xj) F(xj 1) =p(xj)whereF(x0) is defined to be 0. So, the probability thatIj= 1 is same as the probabilitythatX=xj, and this can be used to generate from the distribution ofX. As an example,suppose thatXtakes values inS={1,2,3}with probability mass function defined bythe following table:p(x)xp11p22p33To 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 intervalthe variable falls into.

3 The following R code does this, and checks the results forp1=.4,p2=.25, andp3=.35:# n is the sample size# p is a 3-length vector containing the corresponding probabilitiesrX <- function(n, p){# generate the underlying uniform(0,1) variablesU <- runif(n)# storage for the random variables, XX <- rep(0,n)# which U are in the first intervalw1 <- which(U <= p[1])X[w1] <- 1# which U are in the second intervalw2 <- which( (U > p[1]) & (U < sum(p[1:2])) )X[w2] <- 2# which U are in the third intervalw3 <- which( U > sum(p[1:2]) )X[w3] <- 3return(X)}X = rX(10000, c(.))

4 4, .25, .35))mean(X == 1)[1] (X == 2)[1] (X == 3)[1] empirical probabilities appear to agree with the true values. Any discrete randomvariable with a finite sample space can be generated analogously, although the use of afor loop will be necessary when the number of intervals to check is large. The processwill be similar when the variable has an infinite sample space one example of this is thePoisson distribution. The probability mass function for the poisson with parameter hasthe formp(x) =e xx!whose sample space is all non-negative integers.

5 The following R program generates fromthis distribution and compared the empirical mass function with the true mass functionfor = 4:# function to calculate poission CDF at an integer x# when lambda = L. This does the same thing as ppois() <- function(x, L){# grid of values at which to calculate the mass functionv <- seq(0, x, by=1)# return CDF valuereturn( sum( exp(-L)*(L^v)/factorial(v) ) )}# n: sample size, L: <- function(n, L){U <- runif(n)X <- rep(0,n)# loop through each uniformfor(i in 1:n){# first check if you are in the first intervalif(U[i] < (0,L)){X[i] <- 0} else{# while loop to determine which subinterval,I, you are in# terminated when B = TRUEB = FALSEI = 0while(B == FALSE){# the interval to checkint <- c( (I, L), (I+1,L) )# see if the uniform is in that intervalif( (U[i] > int[1]) & (U[i] < int[2]) ){# if so, quit the while loop and store the valueX[i] <- I+1B = TRUE} else{# If not, continue the while loop and increase I by 1I=I+1}}}}return(X)}# generate 1000 Pois(4) random variablesV = (1000, 4)

6 # empirical mass functionMf <- 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 itb <- 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 functionpoints(b, dpois(b,4), col=4)lines(b, dpois(b,4), col=4)2 Sampling from continuous distributionsContinuous random variables are (informally) those whose sample space is composedof real intervals not exclusively containing integers. In a continuous distribution the prob-ability of taking on any particular value in the sample space in 0; probabilities can only beassigned to intervals in a continuous distribution.

7 For this reason the logic of the previoussection does not apply directly and other methods must be The inversion methodIt is a fact that ifXhas CDFF, thenF(X) has a Uniform(0,1) distribution. The proofof 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))=xSo the CDF ofF(X) isx, which is the same as the CDF of as Uniform(0,1). HereF 1denotes the inverse of the CDF (also called the quantile function) and is defined as thefunction which satisfies:F(x) =y x=F 1(y)In other words,F(F 1(y)) =yandF 1(F(x)) =x.

8 The fact above implies that ifXhasCDFF, thenF 1(U) will have CDFF. So, if you are able to calculateF 1, and cangenerate uniforms, then you can generate a sample an example, suppose you want to generate from a distribution with CDFF(x) =x1 +xforx (0, ). This distribution has mean and variance equal to . To calculateF 1(y),you specify a value foryand solve forx:y=x1 +x(1 +x)y=xy=x xyy=x(1 y)yy 1=xsoF 1(y) =yy 1. ThereforeU1 Uwill have CDFF:# n is the sample <- function(n){# generate uniformsU <- runif(n)# return F^-1(U)return( U/(1-U) )}X <- (1000)# empirical CDFv <- seq(0, 20, length=1000) <- rep(0, 1000)for(i in 1:1000) [i] <- mean(X <= v[i])# true <- v/(1+v)plot(v, , xlab="X", ylab="F(X)", main="Empirical vs True CDF", col=2,type="l")lines(v, ,col=4)Due to the long tails, this distribution is a good candidate for a trial distribution inrejection Sampling , which we will mention later.

9 As a second example supposeXhasCDFF(x) =(11 +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 bitmore involved:y=(11 +e x) 1/y= (1 +e x) (1/y)1/ = 1 +e x log((1/y)1/ 1)=xSoF 1(y) = log((1/y)1/ 1).# generate n samples from the skew logistic# distribution with parameter <- function(n, Theta){U <- runif(n)return( -log( (1/U)^(1/Theta) - 1) )}X <- (1000, 4)# empirical CDFv <- seq(-10, 10, length=1000) <- rep(0, 1000)for(i in 1:1000) [i] <- mean(X <= v[i])# true <- (1 + exp(-v))^(-4)plot(v, , xlab="X", ylab="F(X)", main="Empirical vs True CDF", col=2,type="l")lines(v, ,col=4)In many cases you can not symbolically invert the CDF (the normal distribution is anexample of this).

10 In principle, you can still use this method in such situations, but youwill have to numerically calculate the quantile function. When we get to the section ofthe course dealing with line searching (optimization), we will do an example like Transformation methodIn some situations you can know mathematically that a particular function of a randomvariable has a certain distribution. On homework 1, problem 2 you were given an exampleof this the transformation introduced there is called the Box-Muller transform, and isactually the preferred way to generate normally distributed such relationships are know, it gives a simple way of generating from a this example, suppose we wish to generate from the exponential( ) distribution, andonly have access to a computer which generates numbers from the skew logistic distribu-tion, and do not know the inversion method.