I FEB 1 9 997 Uus. 4As UIRttAR RESEARCH SUPPORT UNIVERSITY OF MICHIGAN BUSINESS SCHOOL DECEMBER 1996 SAMPLING TRUNCATED POISSON AND MULTIVARIATE NORMAL DENSITIES VIA THE GIBBS SAMPLER WORKING PAPER #9612-40 CRAIG CUMIBUS KIMBERLY-CLARK CORPORATION AND PAUL DAMIEN UNIVERSITY OP MICHIGAN BUSINESS SCHOOL AND STEPHEN WALKER IMPERIAL COLLEGE, LONDON, UK

Sampling Truncated Poisson and Multivariate Normal Densities via the Gibbs Sampler Craig Cunlbus Kimberlv-Clark Corporation. Atlanta. GA. USA Paul Damien School of Business. University of Michigan. Ann Arbor. MI, USA Stephen Walker Department of Mathematics, Imperial College, London. UK. ABSTRACT Gelfand et al. (1992) consider the Bayesian analysis of constrained parameter and truncated data problems within a Gibbs sampling framework. In this paper we consider the truncated data problem and concentrate on sampling truncated Poisson and multivariate normal densities within a Gibbs sampler. The proposed method bypasses the need for rejection sampling and/or algorithms such as the Metropolis-Hastings and sampling-resampling. In the truncated multivariate normal case we introduce a Gibbs sampler in which all the conditional distributions can be sampled via uniform variates. whereas the method of Robert (1995), which also runs on a Gibbs sampler, has conditional distributions requiring rejection algorithms to obtain random variates. Key Words: Latent variables, Gibbs sampler, uniform random variables. 1

I. Introduction Let f be a continuous density function defined on the real line. The solution of generating a random variate X from f. starting with the assumption that it is possible to sample uniform random variables from the interval (0. 1) to provide an iid sequence of such variables. was developed in Damien and Walker (1996 a). To motivate their central result stated as a Theorem later. let f(x.) x exp(-x3)l(a < x < b). where > ' < a < b < x and I represents the indicator function: write such a (...ity as ~(3. a.b). Let f(.x) x (1 - x)3-. where 3 > 0 and 0 < a < b < 1: write such a density as B(.3. a. b). Finally, let L'(a. b) denote the uniform distribution on (a. b). It is obvious that both ~(S3,a. b) and B(3. a, ) can be sampled via the inverse transform method. Now, the basic idea is to introduce a latent variable'Y. construct the joint density of Y and X with marginal density for X given by f. and to use the Gibbs sampler (Smith and Roberts. 1993) to generate random variates from f. In particular, Damien and Walker (1996 a) show that all the conditional distributions in such a Gibbs sampler will be one of the three distributions. U(a, 6). E(3, a, b) or B(/3, a, b). To this end they prove: Theorem (Damien and Walker 1996 a). If f(x) x fi i(x), (=1 where the gj are nonnegative invertible functions (not necessarily) densities: that is. if gf(x) > y then it is possible to obtain the set Al(y) = {x: g(x) > y}. then it is possible to implement a Gibbs sampler for generating random variates from f in which all the full conditionals are uniform distributions. Damien and Walker (1996 a) exemplify the use of this Theorem by sampling all univariate continuous densities that appear in Johnson and Kotz. Also they illustrate the method for several Bayesian nonconjugate models and nonparametric models. Damien and Walker (1996 b), Walker and Damien ( 1996 a b) apply the above Theorem to provide a full Bayesian analysis of circular data. neutral to the right processes, and mixtures of Dirichlet processes. respectively. 9

I Since an indicator function satisfies the conditions of a g, it is no more lifficult to sample a truncated or constrained randlom variate using the approach of Dainien and Walker than in the nontruncated or unconstr Wined case. Therefore. if the sampling from f. or a truncated version of f. is required within the context of a Gibbs sampler the Theorem provides a means of doing this via the introduction of latent variables which has the effect of extending the Gibbs loop (i.e. full conditionals) by the number of latent variables introduced. In this paper. the Theorem is used to generate random variates from truncated Poisson (Section II) and multivariate normal densities (Section III). A brief discussion in Section IV concludes the paper. II Truncated Poisson density We consider the random variate generation of a left truncated Poisson variable required within the context of a Gibbs sampler which is 'running' over a sequence of full conditionals. one of which turns out to be a truncated Poisson, for example, in a truncated data problem involving a Poisson model (Gelfand et al.. 1992). Devroye (pg. 489) identifies evaluating factorials as a problem with the sampling of a Poisson variate. Our approach by-passes the need to evaluate factorials. To introduce our approach. let X be a Poisson random variable with fir x AX/.r! I(.x E ) where f, = P(X = x). A > 0, and l = {0. 1.2,...}. Following the Theorem we introduce two latent variables Y and Z such that the joint density with X is given by fx,y,z(x,y,) (y <,\, x< x. The resulting full conditionals are given by YI(X = x) - U{Ox) Zl(X = ) -U (0,!-1) XI(Y = y, Z = ) - U ({n.(y),.... n(y, )}) where, if A > 1 then n(y) is the smallest integer greater than logy/ logA if y > 1, and n(y) = 0 if y < 1, else if A < 1 then n(y) = 0. If A > 1 3

then m(y. ) is the largest integer m.(z) such that m(z)! < 1/-. else if A < 1 then mn(y. ) = rnin{mn(z). a(y)}. where m(y) is the largest integer less than log y/ log A\. Note we can take - = r/x! where r -- '(0. 1) leading to n(z) =int{r-4 (r-lr + 1)) - l}. where r denotes the gamma function. Perhaps a less computational approach to obtaining mn(z) is that mz(z) is the largest integer for which m(z) S log k < - log r, k=x+l and so there is no need for us to evaluate a factorial. Suppose then we need to sample from 1f.x AX/x! I(x E D) where D = {a. a + I...}. The algorithm is now given by the Gibbs sampler with full conditionals Yj( =.) U (0AT) Z (X =_ r) - U' (0.x!-1) X(Y = y, Z = ) ~ U ((max{a, n(y)},..., m(y, )}). and so is essentially no more difficult than the non-truncated case. III. Truncated multivariate normal density In this Section we describe a method for sampling truncated multivariate normal variables using a Gibbs sampler. To introduce our ideas we consider the standard univariate normal distribution. Let X - Ni(0, 1), that is, fx(x) x exp (-x/2). Introduce the latent variable Y which has joint density with X given by fx,,Y(, y).x I(o.,exp(-x2/2))(Y)We then have the following full conditional densities: YJ(X = x) U (0,exp(-x2/2)) XJ(Y = y) (-/-2 log y, -2 log y). 4

Since we are only going to be interested in the. samples we can define I = rexp(-.~2/2) where r - (0 1) and therefore we can construct the sequence {X'.1 } ''A A\1+1\(Xn = -V. rTn) / (- 2 - 2logr /r - B -log) g where the r, are iid '(0. 1) and also independent of the A.,. Obviously. as n1 - x: then ~ — n 4.V1(0. I). Implementing this idea for a truncated standard normal variable is no more complicated. Suppose we wish to sample from the density given by f(.r) x exp (-r/2) I(r E (a.b)). Our required algorithm is now easily seen to be given by the sequence Xn + I (Xn = X. rn) ( T (max{a,-r2 2 log r }, min{b., x2 - log r }). Robert (1995) proposes a rejection based algorithm for sampling a truncated univariate normal distribution. We anticipate that the vast majority of cases involving sampling truncated normal distributions would arise in the context of a Gibbs sampler. so our algorithm is effectively extending each Gibbs loop by only one more full conditional, which happens to be a uniform distribution. Robert (1995) then proposes a Gibbs sampler for sampling from a truncated multivariate normal distribution. We can greatly simplify his algorithm using our latent variable idea. Therefore, consider./.....v,(x,..., z) cc exp (-1/2(x - a)' -( -,u)) I(x E A). where we assume, as does Robert, that the bounds for xi given xi are available and given by, say, (ai, bi). Therefore fxyx(xilx-g) x exp (-i1/2(x - Vi)2/a?) I (xi e (ai, bi)) for i = 1..... p. are the full conditionals. and vi = pi - jAi(j — ili)ij/eii and a'7 = l/ei. where ej is the ijth element of S-1. Robert uses his rejection algorithm for sampling these truncated univariate normal densities. However, since we are already in a Gibbs sampler it seems appropriate to implement 5

the latent variable idea. \Ve do not need to introduce p latent variables. one is.sifficient. \e define the joint density of (.~....... 1) by.v....v..-(. l....,. y) xcexp(-y//2)I ( > (.r - t)l'-1( - t)) I(.r.). The full conditional distributions are given by f', x_.A'(.XlIr-,.y) x I (i G Ai). where Ai = (ai, hi) n B, and B, is the set {xl._-: (x - l)'E-l( - p) < y} and so the bounds for. B are obtained by solving a quadratic equation. The full conditional for (Y jX) is clearly a truncated exponential distribution which can be sampled using the cdf inversion technique. Therefore we have a Gibbs sampler which runs on p + 1 full conditionals which can all be sampled directly using uniform variates replacing the;p full conditionals of Robert which are sampled via rejection algorithms. E.arnple Truncated N2(0. E). We consider the example presented by Robert (1995) which involves sampling a truncated V2(0. E) distribution, where -=(p 1) so v._- (1 _p)- 1 -P(l-P2)-p I - - -p(l-p)- (1 - p2)- and A4 is the ball centered on (71, 2) with radius r. Therefore i = /, - V/r2 -Q(- -x)'2 and bi = 7; + r2- (-y - j)2. and B, - (xjp- - )( - ), x+ (y -. )( l- p2)) 6

~_1 Figure 1: Plot of 10.000 samples from truncated N2(0.S) distribution obtained from Gibbs sampler for i(j) = 1(2).2(1). Note that y > xI since (YI.JliXX2) is an exponential distribution restricted to the set ([(x - 2pxlx2 + x.]/(1 -p2), ) and x?( 1 - p'2) < 1 - 2ipxlx2 + x2, for i = 1,2. We performed the Gibbs sampler for the truncated bivariate normal distribution given and took p = 0.9. (71,22) = (1/2. 1/2), and r = 1. We took 10.000 samples with computing time of 5 seconds and a scatter plot of the samples appears in Figure 1. IV. Extensions and conclusions In this paper we have proposed and illustrated a method for generating random variates from truncated Poisson and multivariate normal densities. The main advantage of the method is that it obviates the difficulties associated with alternative methods such as rejection algorithms (Devroye. 1986). Metropolis-Hastings (Tierney. 1994; Chib and Greenberg, 1996) and sampling-resampling (Smith and Gelfand, 1992). These algorithms require identifying dominating densities, calculating supremums, and acceptance rates: these may be difficult to obtain in many contexts. Required within the context of a Gibbs sampler our algorithms are particularly appropriate since the method leads to the addition of one or two uniform full conditionals 7

I into the original Gibbs loop. \\e restricted ourselves to the Poisson density since simulating froim other discrete dellsities. using our method. is straightforward. Likewise. it is no more difficult to simulate other multivariate densities such as the Student-t. using our approach. As examples, we give here the full conditionals to sample the geometric and negative binomial distributions. 1. Geometric. Here we have f, x (- p)' where p E (0.1) and x E {0. 1...}. The joint density of interest is given by fJy(x.y).x I(y < (I- p)Y). leading to the full conditionals fw'lx(ylx) = C(0, (1 - p)). and fxy-(xIy) = t ({0.1...., n(y)}), where n(y) = int[log y/log(1 - p)]. 2. Negative binomial. Here we have fx oc (x - 1)!(1 - p)I/(x - k)! for p E (0. 1) and x. E {k, k + 1.....The joint density of interest is given by f\x.Y.z.(.r. y.:. w) x l(y < (x - 1)!, =(x - k)! < 1, w < (1 - p)), leading to the full conditonals fwix(wui) = U (0, (1 - p)), fylx(yx) = U (0, (X - 1)!), fzix(Zlx)- = U (o0, (x- k)!-), and fxlwvYv.z(x w, y,) = U({n(y),..., min{n(w), n(z)}}) where n(y) is the smallest integer such that (n(y)!-1) > y, n(z) is the largest integer > k such that (n(-) - k-)! < 1/z, and n(wz) = int[logw/log( - p)]. As with the Poisson case we can circumvent the need to evaluate factorials. 8

References (hib. S. andl Greenberg. E. (1995) Understanding the Metropolis-Hastings algorithm. The. American Statistician. 49, 327-335. Damien P. and Walker S.G. (1996 a) Sampling probability densities via uniform random variables and a Gibbs sampler. Submitted for publication. Damien P.and Walker S.G. (1996 b) A full Bayesian analysis of circular data using von Mises distribution. Submitted for publication. Devroye. L. (1986).Non-uniform random variate generation. Springer. New York. Gelfand.A.E.,Smith.A.FM. and Lee.T.M. (1992). Bayesian analysis of constrained parameter and truncated data problems using Gibbs sampling. Journal of the.American Statistical Association 87, 523-532. Robert. C.P. (199.5). Simulation of truncated normal variables. Statistics and Computing 5, 121-125. Smith. A.F.M[., and Gelfand. A.E. (1992) Bayesian statistics without tears: a sampling-resampling perspective. The American Statistician. 46. 84-88. Smith. A.F.M., and Roberts, G.O. (1993) Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods. Journal of thE Royal Statistical Society. Series B. 55, 3-24. Tierney, L. (1994) Markov chains for exploring posterior distributions (with discussion). The Annals of Statistics. 22. 1701-1762. Walker. S.G. and Damien. P. (1996 a) A full Bayesian nonparametric analysis involving a neutral to the right process. Submitted for publication. 9

I Kresge Business Administration Library Stacks WP5300 1996 no.9612-40 Sampling truncated Poisson and multivariate normal densities via the Gibbs sampler WValker. S.G. and Damien. P. (1996 b) Sampling a Dirichlet process mixtl're modlel. Sub.roitted for publication. UNIVERSITY OF MICHIGAN BUSINESS SCHOOL KRESGE BUSINESS ADMINISTRATION LIBRARY 701 Tappan, K3330 Ann Arbor, MI 48109-1234 DATE DUE SEP 11998 10