I FEB 1 8 1997 'AKESGs U. ASM& LIBRARy RESEARCH SUPPORT UNIVERSITY OF MICHIGAN BUSINESS SCHOOL DECEMBER 1996 SAMPLING NONSTANDARD DISTRIBIUTIONS VIA THE GIBBS SAMPLER WORKING PAPER #9612-34 PAUL DAMIEN UNIVERSITY OF MICHIGAN BUSINESS SCHOOL AND STEPHEN WALKER IMPERIAL COLLEGE, 180 QUEEN'S GATE, LONDON

I Sampling Nonstandard distributions via the Gibbs Sampler Paul Damienl and Stephen Walker2 l)epartment of Statistics andl Management Science. School of Business. University of.\lichligan. Ann Arbor. 4819- 1234. 'S A. D> lepartment of Mat lieiatlics. Imperial College. 180 Queein's (;ate. Lonidon S\\T 2BZ. S.%u1iilmary SiIlmpleand fast atgoritllns to iample i)u' Itandard dlistriluti ions as a pp)earilg in Devroye ( 198i) are d'evelul'ed. [They replace the need( for rejection based appl)roac'iles and are particularly ap)l)ropriate for use in a (hih))bs sampling conIt\f I/ I l.>: Latent varial*|l>. (;ibbs) sarmpler. I litormI (distribittion. I

1 1 Introduction leI / 1 a llt illllotts de it l itv iinctiol l itell'dl on I ' re';al lille'. I t1,' -oiIi iot l, 2eit''Ite t hi it I' i ndo.ll variate. lro'm /. sta.lt'tiIIg it I t <I iL'11111IptI 1l I lia it is possiblek to sample i' tiform randoItll variable.- from ihe inll 'r-lt 1 ) 0. 1It providle a ii l st'(l.l ence of -sl 'ch xvaria lt'. w\as }e'(velop)ed iiI 1)IItiIiti aU d(l \\V Ilker l }I 9 6). 1to Iiloti ale' ll hiir r'( sitI slateil;is a Ilt,o ur i ialatr. l e (.'t t x ('x p -.1' \/I( <.' < l). \W('{r( i > 0. 0 <., < h <. x nli( / r(')irs("ilits I II ll(licator 1't i icltill: r\\'ite('. itcl a, 'i\ i- 3 (. a, 1..h. l, /I.f. \x (1 -.r". whei. e i > () nl 0 K< a < A, 1: writ.e -<clt a,le'ti.ilv it.- i a. h.). It is ohlviotIs tliat tlhese detlsities rat Ine "a-itllt'fl X'ii! It' iInvere t I fraliisori tIti li(lud. Fitiallv. let 1(. h) (ldelt)t- t ltet uitlit 'il rli.sI rillitioll on ((/. ). Now. thle )basic idea is to introdur l a laoitnt 'variall(e. 5. 'tist,;tut tlei joint dellnsity of I and X witlh mlarginal density 'or A.\ X ivei' l}v. an ti to mise lie (;ibs samtpll1er (.Snith iand Roberts. IW99)) It) yieier-;te r'tindoin variates tioin. Il lparticular. Datnien and Walker ( I99i shlow tlin aill te cotndlitional dlistril)btions in a snutt a G;ibbsl samll)ler will 1)' onl.,tf tlie' 1 itee dist rilO>ti(tots. ( 10. b). 1{ 1. i. f). e(.. 0I. ). To( t lis (end li t - iprove: Theorem ( D)aini' and Walker. 1 )96). If L J.',') x J!t(x.). '=1 'v. Icre I lie qg are nontneI ativ e iIve('rt ile fIuntiiotns ( 1iot lnect'ssatril ) ldetsiies -- that it i ts. i,x.) > 1/ tltetl it is possible tot obltaiI t he se-t. It//) = {./: f/.( ) > //} ie(CtI it ils possi)le to i1Il)lti-lieent a libbs sampller for getteratilla g v'ariat(e< i'<tln / int wlicli all bllt one conditional dlisTrilitions will unilform dist riribI)atmien and \\alker I l!)1fi exe'iJlplitfy - h tse of tlhe liTheorem bLv sam)tijilo,.dli llivariat'e t)tliltuoltl > tleniiti'.s thlat appear it.JohIso)t and l\otz. A.lso it [' ill -Ist rat Ilie t it 't t l l t ' t ievera v l 'iil I 11an lon<t'otjnaie dt m lte ls;Iind tiil>araet tlric notde'ls. ('Cttl)i-. )alienl and Walker (I99() illtstlrate thlle -,', t I i ' li lit'orenm to( sitimilalte discre(te andtl tI 'tincat ('led i lt ivariat e dcletsilt ifes..1d(1 \\Walker antld Dam)ieni ( ' )(i) ilse t lie method to sallple ('om a Ilixt irt'( o D)irichle't process model. -

I \\W' c'sisidlter toulr.of tof ' l )tl.otnsi aul(ard (li tril)Iiol.l- s ('oili0.lere;tl iIl I 'lchapli t ') *1' I ),\I'r \'<V' 1!i 1)() i mIlc' l ('ii } l) l i sm11pIlt'(l dir' ctl't. 'I'v l Zi!' li.i ri Fl.*'in i'- 1l'e.'l ilt iingulistiC- ti l.d i'il sc.'i'll. w\ilini' llk'!]il', list il l iul i- *'rIni i errle in t\ m |l >i l -lciersinitre. It is kio\n- tI ha)l I)rcc.l 1)11' (cn! *tii'int |lv sitinilat, a Z/ipf rinii(.ldoIi variale,. il is.triai,'lli',1wt iardI t)i, I~irr;,lF,I l'laic'k vari ' as well. 'Iwto oilier well known, li.-. ' ilt ii-. '<' lin l'tieialistti invx'rs't (i;1 ailhti aI ill. eir.soni \' listril, it io is. In I iis plapler Iltew alorithlts to sai,)ple thle.se (listribll) io)liI anrc',1'I'Iiit I'. I li~'<, algorithlns are. all alteriative to tlie rjiecttion algorillinus whitCl nit lan.ar i; l)Dtv'rot' 1 S198(i. andt i I ba.-ed onl tile 'I'lureni t111 1ivent arli1er. 2 The Zipf distribution ih 1t' ic'INanIill zeta t uI ict1o is, gin\'ti 1)V; 1 ((,) = wltre t - =.s + j/..s. f are real nttnllbers. ald j is tile s<(quare root of litislltS Irn'. Simple expressions fo t lie zeta function are know' n 1. P'ifl rase.s. For e'Ixa ti)leC if,I is art /Ift /f r' lien ( 2i )= - " /l,. (2. )!,vl....v H, is tlbe atth Bernwolli nxlber ' 'Filc'marslh. 1151. p20). I nll Zipl) distrilil) ion las onei parameter (I > I. and is definled y I he prol)a1liiit it's '=!-. ti> 1). FI llIowing tlhe Tlieoremi. tlit (ill sa Ipler is lased on t lie jo iIt dvi stiib.u rtio Jiv I'll 111 to p)rol)ort i)nalitv liy n.:. \ 11() < i'f < I). I l'el mtI ' 1titll r()o tldil i loals are m ix' ' ve, l,uij = I0. i' ").:3

I \t. lit~ RJ = t; ir; 1 1ft ii.;)1 —ji =./fJJ 3=' 1,,'i 'li; F,. \\ '-11 C(I t- ri t l t I.arkov chain A ' 1 il..... 11 t li t 1. V., /;,. — X.v,_,A\.. ]=- t'(I.I2....,,., }, wultr Il,. is railot tlld.aiv(-' 3'! ilalt[.\',r '"] wlic'rt' \l;e.:. t;it'? ii, /1 i(j. tlnl iii(ldepeldel t (' t '.\,,.! re illt [.j (1denotes t1( ' lair('st iIIl(Tr If-. aitll],I ('(1 1t ll to.r. \\' rlanl tle above aloorithm for a 1luniber of valie, lu;) r ( a1 i1i ('I(a' a.M o' (lt'rt(1 20.0()00 raIndom variates. (.onputing tim ti fi'Or ('acli ri ll was ail)lO.xiiltIately 10 seconds anil( the (lstimate.s for t l r(e rcil)roal Wf tl. I't, 1',tiuitt at I 2. I and it are. respectively. 0.0t). 0.923 and 11.1. i.\ \A Mont ('iri iaip)roximilatiOl to Ill i' e lil [a'iictioll is given Ibv,,, = 1 (.A.-i). A-=] d(1 N\. N-....... Iar' lit'.sample]('.- ob)tained filo] the Ziplt' listril titon.. 'lT(h (exa.t V'alles. are (i/;,.!()'^ an 4. 'l/,.' wlhich are ill od aoreelliill with *)I I' t( I ti lla (s. 3 The Planck distribution I lle Planck is ita \\w param-ter dlistrilut ion with deinsitv, givenl bv /'t. = -..+1 [I I/ + I ) (/, + 1 ) (" - I '.h(ere.r > 0. (I > I.is i 'a' p)arameteer. I > (1 is a scale )paraniftier. Our;jipproaich Iin tlihe T orem 1 o \lvest! conisiderat ion of tl lie jit (lelsity giWell 1|> v, I['trl[)rtfioialtit- lp'v.l(..'t.-'x I\/(J <.t <. r". t'(exp(.r) - 1 < I i. I

1 in a whliere f. 4 > 1). ( Clearly liet' mara-ilal dl rib t tiui ol. i>. as r(tlii 'ui. I''li( bi'll * 1ti1itiuiial 'list ril)ltion1 s arn all iiiitiwm ald,uive'i 1x'v./''./J!. ii = -I ' " Ia.. t, I L \\W* tl t li.- al&arit 1il witlf -- - '2 atnd 4 = -.->.and cu ltle'dl 1lf iti.aill-,lt'> fr.tlr tle i Placck dlistriil)itioul. ('OilpI inntg ltill' wat 1:3 *')l-l andl t tist ut.iail reilre.seiitatiolI o(f I lie.saml ple is Ipreseit(e(1 it 1. igiitI 1. 1 - iI-t're 1: flistoogram repre-sentation ofl samples obtlained tfrom Ihe Plla'ck 1di.-K riltion wvith paIrai(eters (/ =2 and i 4= 10.3. 4 The generalised inverse Gaussian distribution O()r alg,1rit ltii i.s pIart icu'larly;il)l)ppropriate 01' t le gelleratlised inverse (;lssiami I (; [( 1 ist il)tit oll. Ac\'cordlit, to De'vrove IIt ere are two dlifficult ies assuoiatt'ei W'ilt I a rejection based algorithI tIo.sanmple the (;GIG distribltionl: firstlv. it r:,.)

I i> uired t c(t i It) alcutlat'e tIle 111(dified i (' s^ etl finiction t tl tl iirdil r.i k;ind: i. -ci',l,1lv. trhei t'Nx)c('tedI I11imhlber F itepal atiols in tI r},eject it)l alturii liII iii liri'e. Tit1' (;I(; di.-tri)lir ioll lIas tiliree parinltlet'lt s \vii II t('d il lll('l iti ti\'lVel IlI) It) i ('i ISt} iit of )'porJti l itot lil l)y /'}.'N r- 'N1 ) ( - -O.-5( /.' -.rji. '.'l'tlf'.x > D(..\ '_ -:. C. and tl.1) > 01. NotI' tli it Ill(',.lillilit llm ll ii\'(' 1t' (;al.iaii (lisl rilii it i s ar ' ai si)(1cial C'as t' oft I l (( i.; \\'(. dcli'i Ilic.,iiitt i, *.il' iil.//t.t.t'. tt' i )' f'.i t. t..,. t. r x I (,< < e(p. " X < 'x(-0. /.'). a' < ('Xl)(-0..t3. ). I 1l t'ltl condlitionals aret iv'tll I)'./'( f t'.. '. t '. = U - '( -O."- I. t "i '"..t. j = '(0). expl( -0(.5/.) 1. /i l.r. t l. t' ) = U '(0. exp(-J0.1. r) ). f - (llm {-/(2 log e)../l } - -/() log ) if A > I./'j... /.i t.. ) = I ( /{(2 log ). -2/) log t) if A 1 t2 (- kg:'( tlo ). mllin{l{p/f lb)l- )2/lo^t}) if' A\ < 1. 5 The Pearson IV distribution I l't*1 'ar..oll I'aiiil o(' (lisi riiuit,it.. has 12 i(indlI)eXrs andr all I)l one can \I, -aMilllI(l, lir'(ctly l si1.Il -t anidlardl dilrillitionls. Tlie excel)tlion is the lPearsoii I\' lisi ril>iltioln..iivten 1l1) u it I'to"i ll i )of propo1) lrtionality })./'(.r) \ (I -r-.f t'x ) 'xp ( - arctan(.r/(t)) we'i-re.r is delftied on -_2c.- -c). > 0. h > 1/2 and c is a real. We defillne I 11'i joilnt deiisit v /.t'. '-.'j) x I(f '< I -- 4.., ) '. - <.?xp)(-t'- arctal i.r/,/))

I Tl' t1,1 l t'onditonals are "ivex ll b).l //'j'..r i / ' (. i -. i../'t.:..r./i-1 / (U. Xl) -p carrcia I. I..1 tillud. n =l ( n(iax v -. \;1l'1' \\l 1a\'-e atssliiiiie'.l \ 11itl -it lo. ofs tetl eraltvt tlitt r > 1). 6 Conclusions I l tli. n [)iaper. w e Il\ave I' i, l i I ew a nd; ii i ll I ay l\\' l i l \ l l t iet atiii iatiiotil variates from t l' Zijl'). Pllic(lk. (;1(. and( Pearsonl\' IV li rilui)t ion-. tli' t hliod.lescril)e(l e(iliiilsi( 111 iore tlhan a (Gibbl)s saiipler involving itniforin till coiditlitional dlistrillltic )iH. olitaiile after the inttrod cl lt iul o latent variAcknowledgements ] Ilis work wva.-s sil)lf>t)i t$i ill part bY a(n EPSR( R()P.\ ".rait. Support for I ravel \was p)roviddl in 1pan l1w ftl' I lit d LKiim(l.-i Engni'l-t lrill." a1Y(1 les'larl'l ( 'tincil and liet I' llix'Vt-i1 of M irli'i.aii. Trili l)pape)r was Coilplelted w\\le. 1i( tii.rt at il or \ \a. \vi-iililtl Imtile'trial ('ollege..Londod n. References ('iii l s... Jaiii. '.. aIIien'.. alker. S.(;. (l!)!)(i) ISaipliliun rincated IPois-tki and liiitivaniatle iurnial eii-til ie -.ii t liet (Gibbs 'alll)l(br. Subittifndf frt!fftl/li'f 1ffl n. D);ainie. P. and \\allkri'.:.(;. * l!'',, iSai|)l)n. p)lrobablit tlilitv,l-',ies via lnIiIIII rIandl(llom \varialbit" dl.1 ( iblr. 'IlInple'r. >'ti/bifitlI Jotr phbfcaf/fion. l)t I'ruv t. [. tl!Ni i \ -.\ a, t-'/,.tf /i{flot, \){ i(f (/; lSt n'l'tion. Spring.er\V't';ia,..New\ ^'trk.

I Kr eso' g BuJs in&ets Adrmini tYatlt.icr Library Ste acks WP,30'A 19i9 no.9612-34 S.alp I irin norn standa rd dl st.r bu. io ti vi th t e.Gibbs< t E1 [r i; F tl tL ~S tizil til. \I.i... 1 1l.. l tIiI*. ( I 1 3!)!9.: ) a13c <'i',lia t i't fijll otil, \l li,' (ij lia -;IIIl)|lIr;,II1l Il.lat. l \[ rkov <'lilii,.\MIl t,' ('<Irlt> 1Icl('ti, 1, \\'ii 1 li.- i ',i s>i!t > 1../,,,B', /f 0/-' ft h h f,{,/!.%' / U/;.-i/,',i.%' /. /t 1//..%', - f.:.. *...).;> -' 1. 'l'li liilarsli. I L.('. l1t) 51 T 791 ' I /' / t, l{!ff /, zf o t/h l"iem-iXf/e1. ().xlirfl I 'li 'ersilv Pr<i^. Oxitord. 19!)1. W\\alker'. S.a.,iuu Daien. i ). I I!))fi) Sa lnpliiig.a Diriihlet l),(',S N1 ixt rltr,,ii'!l. f'inl S.;:~, 8:i,-,lotr l tblicatiotl.