In this short paper we summarize the computational steps of Adaptive Radial-Based  Direction Sampling (ARDS), which can be used for Bayesian analysis of ill behaved target densities. We consider one simulation experiment in order to illustrate the good performance of ARDS relative to the independence chain MH algorithm and importance sampling

Bauwens, L. (Luc)

Bos, C.S. (Charles)

Dijk, H.K. (Herman) van

Oest, R.D. (Rutger) van

Bauwens, L

Bos, CS (Charles)

van Dijk, Herman

van Oest, RD (Rutger)

EUR Research Repository

Explaining adaptive radial-based direction sampling

Bos, C.S.

Bauwens, L.J.

van Dijk, H.K.

van Oest, R.D.

VU Research Portal

Explaining Adaptive Radial-based Direction Sampling

NARCIS 

Erasmus University Digital Repository

EXPLAINING ADAPTIVE RADIAL-BASED DIRECTION SAMPLINGEconometric Institute Report EI 2003-37Luc Bauwens, Charles S. Bos, Herman K. van Dijk, and Rutger D. van OestCORE and Department of Economics, Universite´ catholique de LouvainTinbergen Institute and D. of Econometrics & O.R., Free University AmsterdamEconometric and Tinbergen Institutes, Erasmus University RotterdamEconometric and Tinbergen Institutes, Erasmus University RotterdamKeywords: Markov Chain Monte Carlo,importance sampling, radial coordinatesIntroductionIn recent decades Monte Carlo methods, suchas Metropolis-Hastings [MH] (Metropolis etal. 1953, and Hastings 1970), Gibbs sampling(Geman and Geman 1984) and importancesampling, have been applied extensively andsuccessfully within Bayesian analyses of sta-tistical and econometric models. However,although Monte Carlo methods revolutionizedthe applicability of Bayesian inference, theliterature indicates a substantial variation intheir convergence behavior. For instance, theperformance of the Gibbs sampler may behampered by strong correlation in the target (orposterior) distribution. A multimodal targetdensity may pose problems to all methods. Acommon difficulty encountered in all samplersis the choice of an importance or candidatedensity when little is known about the shape ofthe target density.In order to sample effectively from ill behavedtarget densities, we have proposed the classof adaptive radial-based direction sampling(ARDS) methods, see Bauwens et al. (2003).The advantages of the ARDS algorithms arethreefold. Firstly, the algorithms are parsimo-nious in their use of information on the shapeof the target density. Only location and scaleneed to be specified as initial values. Secondly,the algorithms are flexible and robust, that is,they can deal with a large variety of features oftarget distributions, in particular multimodal-ity, strong correlation, extreme skewness andheavy tails. Thirdly, the algorithms can dealwith multiple linear inequality conditions onthe parameter space without any additionalcomplications for the implementation.The ARDS algorithms can be characterizedas ‘black-box’ samplers (Chen et al. 2000).They extend earlier methods like the method ofBox and Muller (1958), the adaptive directionsampling algorithms proposed by Gilks et al.(1994), the mixed integration method by VanDijk et al. (1985), and the spherical integrationmethod by Monahan and Genz (1997).In this paper we briefly outline the com-putational steps of ARDS. For an in-depthdiscussion, the reader is referred to Bauwens etal. (2003).Overview of ARDSMost simulation algorithms for posterior dis-tributions generate random drawings in theoriginal parameter space. However, several re-searchers advocate to simulate in a transformedspace, where the simulation is more efficient insome sense. For example, if there exists strongcorrelation between two random variables, anorthogonalizing transformation reduces serialdependence in a Gibbs sampling scheme. TheARDS algorithms rely on this general idea.They are based on a composite transformationto radial coordinates. For expository purpose,we treat this transformation explicitly in twosteps. The first step is an m-dimensionallocation-scale transformation from x to y,which aims at standardizing the candidatedensity with respect to the location, scale andcorrelations of the target density. The secondstep is a radial transformation, which maps anm-dimensional point y to an (m−1)-dimensionaldirection η and a one-dimensional distance ρ.After applying this composite transformationto a candidate drawing x from the normaldistribution (or any other elliptically-contoureddistribution), the direction η is evaluated byeither an MH step or an importance samplingstep. Next, a new distance ρ is generated,conditional on the direction η, from the exacttarget density by means of the numerical inversetransformation method. Finally, the resultingdirection and distance are transformed to theoriginal space by inverting the composite radialtransformation. This procedure is repeated inmany iterations.The ARDS algorithms have a natural inter-pretation. Basically, in one iteration, first adirection η, which uniquely determines a linein the m-dimensional original parameter space,is generated from an elliptically-contouredcandidate density. This line is either acceptedor rejected in an MH step, or it is assigned animportance weight. Next, a distance ρ, whichuniquely determines the position of one pointon the line, is directly generated from the targetdistribution. As sampling on any given linemimics exactly the target density, ARDS shouldresult in robust Monte Carlo estimates. Werefer to the MH variant of ARDS as adaptiveradial-based Metropolis-Hastings sampling(ARMHS), while we refer to the importancesampling variant as adaptive radial-basedimportance sampling (ARIS).In ARDS, an adaptive procedure is con-sidered to generate more efficient directions,with higher acceptance rates or a more uniformweight distribution. In this adaptive procedure,the location and covariance matrix of thecandidate density are updated by replacingthem by their Monte Carlo estimates. TheseMonte Carlo estimates are obtained from aprevious run of the algorithm. So, ARDS isadaptive in the sense that it aims to improvethe approximation of the candidate distributionto the target distribution by considering severalsampling rounds.The radial transformationThe radial transformation is at the heart ofthe ARDS algorithms. In a general form, theradial transformation from y to (ρ, η) may becharacterized by its inverse transformationy1 = ρ h1(η1, . . . , ηm−1),... (1)ym = ρ hm(η1, . . . , ηm−1),where h(η) = (h1(η), . . . , hm(η)) are differ-entiable functions, see for example Muirhead(1982, Section 1.5). The radial transformationmaps the Cartesian coordinates y to a positionh(η) on the unit circle and a stretching factor ρwhich further determines the position of y givenh(η). It turns out that ARDS can be applied toany transformation satisfying (1).A well-known special case of the generaltransformation described above is the two-dimensional polar transformation from y ∈ R2−1 1ρy2y1ηFigure 1: The relationship between Cartesiancoordinates and polar coordinates in the two-dimensional case.to (ρ, η) ∈ R+ × (0, 2pi). Figure 1 illustratesthe relationship between Cartesian coordinatesand polar coordinates. A feature of the polartransformation is that ρ ∈ R+, implying thatthe direction η only defines one half of a line,and not an entire line. However, in ARDSthe distances ρ are sampled from the targetdensity on entire lines, so that information fromthe whole target density is considered given adirection η. Two other drawbacks of using thepolar transformation in ARDS are the following.Firstly, it is possible but not straightforward togeneralize it to m > 2 dimensions. Secondly,the transformation is computationally not veryefficient. We therefore propose a transformationwhich satisfies (1), is easy to generalize to morethan two dimensions, is more efficient thanthe polar transformation, and which is suchthat the direction η does define an entire line.The proposed transformation maps y to m − 1Cartesian coordinates on the unit circle and astretching factor ρ. This is illustrated in Figure2 for the two-dimensional case. More details aregiven in Bauwens et al. (2003, Subsection 2.1).One iteration of ARDSAn iteration of ARDS consists of 5 steps.Step 1 of the i-th iteration: A candidatedrawing x∗iis generated from a normal candi-date distribution with mean µ and covariancematrix Σ.Step 2: This drawing x∗iis transformedto a candidate direction η∗iand a candidatedistance ρ∗iby applying the composite radialtransformation.−1 1ρy2y1ηρ> 0ρ< 0Figure 2: The relationship between Cartesiancoordinates and radial coordinates in the two-dimensional case.Step 3: An MH step or an importance sam-pling is applied to η∗i, resulting in the directionηi. In ARMHS, this direction is given byηi ={η∗iwith probability α(ηi−1, η∗i)ηi−1 with probability 1− α(ηi−1, η∗i)(2)for some acceptance probability α(ηi−1, η∗i),which is defined below. In ARIS, this directionηi = η∗iis assigned importance weight w(ηi),which is also defined below.Step 4: A distance ρi is sampled from p(ρ|ηi),which is the target density of ρ conditional onthe generated direction ηi.Step 5: The resulting direction ηi and distanceρi are transformed to the original parameterspace by inverting the composite radial trans-formation. This yields one effective realizationxi from the target density.It can be shown that the acceptance probabil-ity in the MH step of ARMHS is given byα(ηi−1, η∗i) = min{w(η∗i)w(ηi−1), 1}, (3)where the importance function w(η), which isalso considered in ARIS, is defined byw(η) =∫∞−∞κ(ρ|η) dρ, (4)and where κ(ρ|η) is a kernel of the conditionaltarget density p(ρ|η), which is defined byp(ρ|η) ∝ κ(ρ|η) = p(x(ρ, η|µ,Σ)) |ρm−1|. (5)Here x(ρ, η|µ,Σ) results from inverting thecomposite radial transformation for the argu-ments ρ and η.In order to obtain the direction ηi, both inARMHS and ARIS, the importance functionw(η) has to be evaluated. Its value can becomputed using a deterministic integration rule,such as the adaptive Simpson’s rule. As thedensity of ρ conditional on η is proportionalto the integrand of w(η), evaluations of thisintegrand, gathered during the deterministicintegration phase, can be used to construct agrid for p(ρ|η). So, using the numerical inversetransformation method, sampling the distanceρ conditional on the direction η is straightfor-ward. We note that one can capitalize on theconstructed grid for p(ρ|η) by generating severaldrawings of ρ for each drawing of η.The integral w(η) has infinite integrationbounds. However, in practice finite integrationbounds have to be considered for its numericalevaluation. In order to obtain these finitebounds, we propose that minimum and maxi-mum values are imposed for each parameter inthe original parameter space. It is often possibleto find sensible bounds, such that practicallyall density mass of ρ given η is included, byeither theory and/or common sense. We notethat our method is in particular suited to dealwith multiple linear inequality restrictions,which sometimes have to be imposed on theparameters. These linear restrictions reduce thefeasible region and hence the integration intervalof w(η), while they do not imply any additionalcomplications for the implementation of ARDS.So, linear restrictions on the parameter spacedo not put a burden on the algorithm, but theymight result in an efficiency gain.An illustrationIn order to illustrate the performance of ARDS,we consider an 8-dimensional trimodal mixturedistribution featuring skewness, some highcorrelations (varying between -0.95 and 0.90)and multimodality. It is given byp1 N(µ1,Σ1) + p2 N(µ2,Σ2) + p3 N(µ3,Σ3),whereµ1 = 1.5 (1, 2, 3, 4, 5, 6, 7, 8)′,µ2 = 1.5 (5, 6, 7, 8, 1, 2, 3, 4)′,µ3 = 1.5 (8, 7, 6, 5, 4, 3, 2, 1)′,Σ1 = Σ2 = Σ3 = I8,p1 = p2 = p3 = 1/3.The marginal densities1 for the eight compo-nents are displayed in Figure 3. It is seen1These marginal densities are constructed using250000 directly sampled drawings.that the target distribution is ill behaved. Weestimate the first and second moments for thisdistribution using ARMHS, ARIS, MH andimportance sampling. This is done in severalsampling rounds. In our adaptive approach,additional sampling rounds are considered aslong as the Mahalanobis distance is larger than0.02. However, in the experiment we allow forat most 8 rounds. In any round, ARMHS andARIS draw 5000 directions and 5 distances perdirection, resulting in a sample of size 25000. Inorder to make the computing times comparable,the MH and importance sampling algorithmsare allowed to collect a larger sample of size250000. The mean for the initial candidatedensity is set at 6 for all eight components.Furthermore, the scale is taken sufficientlylarge so that MH and importance sampling caninitially cover the whole target density.The estimates of the location and scaling pa-rameters are reported in Table 1. Furthermore,the table also contains the “true” parametervalues, obtained from 250000 directly sampleddrawings. It is seen that ARMHS and ARIS doa very good job, whereas MH and importancesampling fail. ARIS only needs 5 rounds toreach convergence (according to our definition),whereas the other three algorithms need themaximum number of rounds. However, aftereight rounds also ARMHS has converged,whereas the other two algorithms clearly havenot, see the reported Mahalanobis distanceswhich concern the final sampling round. Wenote that the average computing times persampling round for the four algorithms arecomparable.The acceptance rates for ARMHS and MH(reported for the final round) show a largedifference in values. Furthermore, it is seenthat the moment estimates, obtained fromimportance sampling, are almost completelydetermined (99.7%) by only 5% of the drawings.In contrast, in ARIS the 5% most influentialdrawings only have 29.2% of the total weight.Finally, in order to see whether the favorableresults are merely a coincidence, we repeat theexperiment above 10 times for different seeds ofthe random number generator. The results arerobust: ARMHS succeeds 9 times and ARISeven succeeds for all 10 repetitions, whereasMH and importance sampling only have successrates of 30%, making the outcome of the lattertwo methods very unreliable.Contact informationCorresponding author: Rutger van Oest, Tin-bergen Institute, Erasmus University Rotter-dam, P.O. Box 1738, NL-3000 DR Rotterdam,The Netherlands (e-mail vanoest@few.eur.nl,tel. +31-10-4088946, FAX +31-10-4089162).The working paper is available as report EI2003-22 athttp://www.few.eur.nl/few/research/pubs/ei/2003/reports.htmAcknowledgementThe authors are indebted to Michel Lubranofor very helpful comments.ReferencesBauwens L., C.S. Bos, H.K. van Dijk andR.D. van Oest (2003), “Adaptive Radial-BasedDirection Sampling: Some Flexible and RobustMonte Carlo Integration Methods”, forthcomingin Journal of Econometrics.Box G. and M. Muller (1958), “A Note on theGeneration of Random Normal Deviates”, An-nals of Mathematical Statistics, 29, 610−611.Chen M.-H., Q.-M. Shao and J.G. Ibrahim(2000), Monte Carlo Methods in Bayesian Com-putation, Springer-Verlag, New York.Geman S. and D. Geman (1984), “Stochastic Re-laxation, Gibbs Distributions, and the BayesianRestoration of Images”, IEEE Transactions onPattern Analysis and Machine Intelligence, 6,721−741.Gilks W., G. Roberts and E. George (1994),“Adaptive Direction Sampling”, The Statisti-cian, 43, 179−189.Hastings W.K. (1970), “Monte Carlo SamplingMethods Using Markov Chains and Their Ap-plications”, “Biometrika”, 57, 97−109.Metropolis N. A.W. Rosenbluth, M.N. Rosen-bluth, A.H. Teller and E. Teller (1953), “Equa-tion of State Calculations by Fast ComputingMachines”, The Journal of Chemical Physics,21, 1087−1092.Monahan J. and A. Genz (1997), “Spherical-Radial Integration Rules for Bayesian Compu-tation”, Journal of the American Statistical As-sociation, 92, 664−674.Muirhead R. (1982), Aspects of MultivariateStatistical Theory, Wiley, New York.Van Dijk H.K., T. Kloek and C.G.E. Boender(1985), “Posterior Moments Computed byMixed Integration”, Journal of Econometrics,29, 3−18.-5 0 5 10 15.1x_1-5 0 5 10 15.1.2 x_2-5 0 5 10 15.1.2 x_3-5 0 5 10 15.1.2 x_4-5 0 5 10 15.1.2 x_5-5 0 5 10 15.1.2 x_6-5 0 5 10 15.1.2 x_7-5 0 5 10 15.1x_8Figure 3: Marginal densities for the trimodal mixture distribution.Table 1: Sampling results for the trimodal mixture distribution.bounds ARDS initialization ARMHS ARIS MH IS truemin. max. mean s.d. mean s.d. mean s.d. mean s.d. mean s.d. mean s.d.x1 -20.00 20.00 6.00 10.00 7.03 4.41 7.07 4.40 1.83 0.04 7.44 0.11 7.00 4.42x2 -20.00 20.00 6.00 10.00 7.54 3.38 7.52 3.43 2.57 0.04 8.10 0.15 7.50 3.40x3 -20.00 20.00 6.00 10.00 8.03 2.73 7.98 2.70 4.09 0.35 10.21 0.08 8.00 2.74x4 -20.00 20.00 6.00 10.00 8.54 2.70 8.47 2.73 6.01 0.39 11.21 0.40 8.50 2.74x5 -20.00 20.00 6.00 10.00 4.99 2.71 5.06 2.70 6.13 0.45 2.96 0.91 5.00 2.74x6 -20.00 20.00 6.00 10.00 5.44 2.77 5.49 2.72 8.91 0.62 2.14 0.61 5.50 2.74x7 -20.00 20.00 6.00 10.00 6.07 3.38 5.96 3.41 10.21 0.32 4.24 0.48 6.00 3.40x8 -20.00 20.00 6.00 10.00 6.48 4.39 6.44 4.46 12.01 1.18 5.14 0.85 6.50 4.42drawings per iteration (η × ρ) 5000×5 5000×5 250000 250000number of iterations 8 5 8 8average time per iteration (in s) 112.9 82.0 110.4 93.9Mahalanobis distance 0.01 0.01 1.57 10.68acceptance rate (in %) 27.8 0.15% most influential weights (in %) 29.2 99.7

Explaining Adaptive Radial-Based Direction Sampling

http://hdl.handle.net/1871.1/8dfc274e-0451-4ce5-8f87-ff6d2c86d6af

Explaining Adaptive Radial-Based Direction Sampling

Abstract

Similar works

Full text

Available Versions

EUR Research Repository

VU Research Portal

EUR Research Repository

NARCIS

Erasmus University Digital Repository