Perfect Gibbs sampling is a method to turn Markov Chain Monte Carlo (MCMC) samplers into exact generators for independent random vectors. We show that a perfect Gibbs sampling algorithm suggested in the literature is not always generating from the correct distribution. (author's abstract)Series: Research Report Series / Department of Statistics and Mathematic

Hörmann, Wolfgang

Leydold, Josef

English

Elektronische Publikationen der Wirtschaftsuniversität Wien

ePubWU Institutional RepositoryWolfgang Hörmann and Josef LeydoldA Note on the Folding CouplerWorking PaperOriginal Citation:Hörmann, Wolfgang and Leydold, Josef (2006) A Note on the Folding Coupler. Research ReportSeries / Department of Statistics and Mathematics, 36. Department of Statistics and Mathematics,WU Vienna University of Economics and Business, Vienna.This version is available at: http://epub.wu.ac.at/1528/Available in ePubWU: May 2006ePubWU, the institutional repository of the WU Vienna University of Economics and Business, isprovided by the University Library and the IT-Services. The aim is to enable open access to thescholarly output of the WU.http://epub.wu.ac.at/A Note on the Folding CouplerWolfgang Ho¨rmann, Josef LeydoldDepartment of Statistics and MathematicsWirtschaftsuniversita¨t WienResearch Report SeriesReport 36May 2006http://statmath.wu-wien.ac.at/A Note on the Folding Coupler ∗Wolfgang Ho¨rmann (hormannw@boun.edu.tr)Department of Statistics and Mathematics, Vienna University of Economics andBusiness Administration, Augasse 2–6, A-1090 Vienna, Austria, andDepartment of Industrial Engineering, Bogazici University Istanbul, 80810Bebek-Istanbul, TurkeyJosef Leydold (leydold@statistik.wu-wien.ac.at)Department of Statistics and Mathematics, Vienna University of Economics andBusiness Administration, Augasse 2–6, A-1090 Vienna, AustriaMay 11, 2006Abstract. Perfect Gibbs sampling is a method to turn Markov Chain Monte Carlo(MCMC) samplers into exact generators for independent random vectors. We showthat a perfect Gibbs sampling algorithm suggested in the literature is not alwaysgenerating from the correct distribution.Keywords: Markov chain Monte Carlo method, perfect Gibbs sampling, couplingfrom the past1. IntroductionMarkov Chain Monte Carlo (MCMC) samplers are very powerful meth-ods for drawing random samples from quite arbitrary distributions. Inparticular they are used in the case of simulations that invoke highdimensional integrals. However, as they produce dependent randomvariables (vectors) they require a convergence assessment. To over-come this problem Propp and Wilson (1996) suggested so called perfectsampling algorithms that allow to decide exactly, when convergence isreached. Although first developed for discrete state spaces perfect sam-pling also can be applied to Markov chains with state space Rd albeitthis is not easy, see (Green and Murdoch, 2000; Wilson, 2000; Murdoch,2000; Murdoch and Meng, 2001).Coupling from the past (CFTP) suggested by Propp and Wilson(1996) is probably the most popular of these perfect sampling algo-rithms. The main building block of all CFTP algorithms is the ran-domizing operation. It is a deterministic function φ taking as input thestate Xt of the chain X at time t and some intrinsic randomness Ut.The randomizing operation returns the new output stateXt+1 = φ(Xt, Ut) .∗ This work was supported by the Austrian Science Foundation (FWF), projectno. P16767-N12.2 W. Ho¨rmann and J. LeydoldRoughly spoken CFTP starts Markov chains from all possible pointsof the state space at some time −T in the past. Using the randomizingoperation with the same intrinsic randomness these chains develop inparallel or coalesce. If at time 0 all have coalesced in a single state thisstate is returned. Otherwise the chains are restarted at some earliertime −T ′ < −T ; see e.g. Wilson (2000) for a short tutorial.Corcoran and Schneider (2003) suggest a new randomizing operationand call it folding coupler. They use it to realize simple perfect Gibbssampling algorithms for different examples. However, it turns out thatthese algorithms are not always sampling from the desired distributionas the claimed monotonicity property of the folding coupler is notfulfilled.The paper is organized as following: Section 2 describes the foldingcoupler and its monotonicity property. In Section 3 the folding coupleris used to compile a Gibbs sampler and it is shown that there aresmall deviations from the target distribution. Section 4 discusses thedifficulties to detect this problem in the design of that sampler.2. The Folding Coupler and its MonotonicityThe folding coupler can be used as a building block for perfect samplingalgorithms. Assume we want to generate numbers uniformly in theinterval (c, d) that is contained in (a, b). Then for a step of the Markovchain the folding coupler draws a trial U uniformly distributed in (a, b).It is accepted if it is falls into the desired interval (c, d), otherwise U is“folded” into the desired interval from the left and from the right. Indetail the folding coupler is defined by the function:f(u, a, b, c, d) =u−ac−a+b−d(d− c) + c for u ∈ (a, c)u for u ∈ (c, d)d− b−uc−a+b−d(d− c) for u ∈ (d, b)The folding coupler is a generalization of the multishift and the multi-scale couplers introduced by Wilson (2000). Its advantage lies in the factthat it can cope with the multiscale-multishift case; one disadvantageof it is that it does not lead to maximal coalescence probabilities, andthus chains need longer time till coalescence than necessary.Much more important is the fact that the folding coupler is notalways monotone in its last two arguments. This is of greatest practi-cal importance as the monotonicity of a coupler allows that only themaximal and the minimal chain must be generated and stored in aperfect sampling algorithm. The problem can be seen quite easily. ForA Note on the Folding Coupler 32 4 6 8 1012345672 4 6 8 101234567Figure 1. Folding couplers for d = 5 and d = 7 as a function of u ∈ [0, 10] (l.h.s.);and folding coupler for a fixed value u = 7.5 as a function of d ∈ [0, 10] (r.h.s.).For both sides a = 0, b = 10, and c = 0.example, fix (a, b) = (0, 10) and c = 0. Then f(u, 0, 10, 0, 7) is largerthan f(u, 0, 10, 0, 5) for values of u between 7 and 8. For these valuesno direct acceptance is possible and thus they are transformed intovalues close to 0. Figure 1 (l.h.s.) shows the graphs of these two foldingcouplers for all possible values u between 0 and 10. It already indicatesthat the folding coupler is not always monotone in its last argument.This becomes even clearer when we look at the r.h.s. of Figure 1. Itshows the result of the folding coupler with a = 0, b = 10 and c = 0for a fixed value of u = 7.5 and d taking all values between 0 and 10.Clearly this function and thus the folding coupler is not monotone ind.Remark. There seems to be no simple way to “repair” this problemof the folding coupler. As we need the uniform distribution there mustbe values that are transformed into 0 if u is larger than c and thisdestroys the monotonicity.It must also be noted that there are parameter that lead to a mono-tone folding coupler. Figure 2 shows such an example. For a = 0 andb = 10 the folding coupler for (c, d) = (2, 7) is below or equal to thefolding coupler for (c, d) = (2, 8). Nevertheless, it is obvious from theFigure 1 that the general monotonicity of the folding coupler in its lasttwo arguments is not fulfilled.3. Using the Folding Coupler for Perfect Gibbs SamplingThe folding coupler samples from the uniform distribution. It can there-fore easily be used to obtain a perfect Gibbs sampler for distributionswith uniform full conditional distributions. As example let us considerthe two-dimensional uniform distribution on the area between density4 W. Ho¨rmann and J. Leydold2 4 6 8 10345678Figure 2. Comparing the folding couplers for d = 7 and d = 8 for a = 0, b = 10,c = 2 and u = (0, 10).and the x1-axis for the standard exponential distribution cut off at 20.That area is is the setS = {(x1, x2)|0 ≤ x1 ≤ 20 and 0 ≤ x2 ≤ exp(−x1)}.The Gibbs sampler for this uniform distribution and a given startingvalue 0 ≤ x2 ≤ 1 draws a random variate X1 uniformly distributed on(0,min(20,− log(x2)), X1 ∼ U(0,min(20,− log(x2)), and then samplesX2 ∼ U(0, exp(−X1)). The next steps repeat this pattern and updatethe previous values of X1 and X2.To turn this very simple Gibbs sampler into a perfect samplingalgorithm it is in theory necessary to run all chains (from all possiblestarting values) in parallel, i.e. using the same randomness for all ofthem. Then the CFTP algorithms start farer and farer in the past ofthe chains till all chains have coalesced to a single value at time 0.Therefore CFTP algorithms need a coupler that allows for coalescence.In other words it must be possible to obtain – when using the same pieceof randomness – the same uniform random number for two uniformU(α, β) distributions with different parameter β. This is the case forthe folding coupler defined above but not the case if we just generatethe uniform distribution with the standard approach.The density of the exponential distribution is monotone decreasing.Thus the Gibbs sampler has the property, that large values of x1 leadto small values of X2 and large values of x2 lead to small values of X1(and vice versa). This property is useful if we can find a coupler thatis monotone for the upper border of the interval. (The lower border isalways 0 so no monotonicity is required for the lower border for thisspecial example.) Then, for a fixed piece of randomness, the generatedvariate of U(0, a) is smaller or equal to the variate of U(0, b) as long asa ≤ b. Such a monotonicity is important as we know that all possiblechains (with arbitrary starting value x2) will always be between thechains starting with x2 = 0 and starting with x2 = 1. It is thereforeA Note on the Folding Coupler 50.2 0.4 0.6 0.8 1-0.1-0.08-0.06-0.04-0.020.020.042.5 5 7.5 10 12.5 15 17.5 20-0.1-0.08-0.06-0.04-0.020.020.04Figure 3. Relative difference between the values generated by the folding couplerand the exact values for the marginal distribution of X1enough to run only these two chains as we know that coalescence of allchains is reached when these two chains have coalesced.Unfortunately the folding coupler is monotone for a large subset ofpossible parameters, but not for all, as we have seen in Sect. 2. Whathappens if we use it for Gibbs sampling as if it were a monotone coupler?We tried this for the uniform distribution over the set S defined above.We ran a CFTP algorithm using the folding coupler and made a chi-square test for the marginal distribution of X1 which should followa standard exponential distribution truncated at 20. The chi-squaretest indicates clear problems for n = 106 but no problems for n =104. Figure 3 displays the relative difference between the histogram ofthe marginal distribution of x1 coordinate of the points generated bythe Gibbs sampler using the folding coupler and the exact (truncatedexponential) distribution. As expected the differences are not large butit clearly indicates that the algorithm does not produce enough valuesvery close to 0 and above 1.The fact that there are only small deviations from the correct dis-tribution needs some explanations. The reason is simply that the algo-rithm is a correct MCMC algorithm converging to the given distribu-tion. It is the monotonicity assumption that is not satisfied and thusthe “perfect” algorithm sometimes returns wrong values as coalescencewas only achieved for the two observed chains, but not for all possiblechains. For illustration purposes we ran the naive Gibbs sampler on Swith a fixed number of 16 updates and with both a uniform startingdistribution and a fixed starting point x1 = 10. For this approximate al-gorithm the generated distribution is so close to the desired distributionthat the chi-square test did not indicate any problems for n = 106.4. Application to the Autoexponential DistributionCorcoran and Schneider (2003) explain how the folding coupler can beused to implement the Gibbs sampler for non-uniform full conditional6 W. Ho¨rmann and J. Leydolddistributions. Following Wilson (2000) they use the slice sampling ideafor this task.To check if the problem with the folding coupler is also visible fortheir example we coded the “folding backward coupling algorithm” onp. 285 of (Corcoran and Schneider, 2003). We only changed the coa-lescence condition as coalescence in one variable automatically impliesthat the second variable coalesces in the next step for a two variableGibbs sampler. We tested the resulting algorithm – as in the paper– for the two dimensional autoexponential distribution with densityproportional topi(x1, x2) = exp(−β1 x1 − β2 x2 − β12 x1 x2)with the parameter values β1 = 2, β2 = 3 and β12 = −1.For the marginal distribution of the first variable we conducted chi-square tests for sample sizes n = 106, 107 and 108 and the number ofclasses close to√n. We used a high-quality multiple-recursive uniformrandom number generator. In this setting problems were only clearlyvisible for n = 108.We can therefore conclude that the difference between the generatedand the exact distribution is so small for this example that the extensiveempirical tests in (Corcoran and Schneider, 2003) did not indicate anyproblem. Thus the fact that the folding coupler is not monotone wasoverlooked.To control our test set-up we applied the same tests to autoexponen-tial random variates generated by naive two-dimensional rejection froma constant hat. The resulting chi-square values were in the expectedinterval even for sample size 108.5. ConclusionWe have shown that the folding coupler does not always fulfill themonotonicity property that is required for the CFTP algorithm pro-posed by Corcoran and Schneider (2003). We demonstrated that thisproblem leads to small deviations from the target distribution. Al-though these deviations are sometimes very small we do not thinkthat the folding coupler algorithm is suitable for practical purposesas it is not clear how close generated distributions are to the desireddistributions for new applications.A Note on the Folding Coupler 7ReferencesCorcoran, J. N. and U. Schneider: 2003, ‘Shift and scale coupling methods for perfectsimulation’. Probability in the Engineering and Informational Sciences 17, 277–303.Green, P. J. and D. J. Murdoch: 2000, ‘Exact sampling for Bayesian inference:towards general purpose algorithms (with discussion)’. In: J. M. Bernardo, J. O.Berger, A. P. Dawid, and A. F. M. Smith (eds.): Monte Carlo Methods, Vol. 6of Bayesian Statistics. Oxford, pp. 301–321, Oxford University Press.Murdoch, D. J.: 2000, ‘Exact sampling for Bayesian inference: unbounded statespace’. In: N. Madras (ed.): Monte Carlo Methods, Vol. 26 of Fields InstituteCommunications. pp. 111–121, Americal Mathematical Society.Murdoch, D. J. and X.-L. Meng: 2001, ‘Towards perfect sampling for Bayesianmixture priors’. In: ISBA 2000, Proceedings. p. ??, ISBA and Eurostat.Propp, J. G. and D. B. Wilson: 1996, ‘Exact sampling with coupled Markov chainsand applications to statistical mechanics’. Random Structures and Algorithms9, 223–252.Wilson, D. B.: 2000, ‘Layered multishift coupling for use in perfect sampling algo-rithms (with a primer on CFTP)’. In: N. Madras (ed.): Monte Carlo Methods,Vol. 26 of Fields Institute Communications. pp. 141–176, American MathematicalSociety.

A Note on the Folding Coupler

Abstract

Similar works

Full text

Available Versions

Elektronische Publikationen der Wirtschaftsuniversität Wien