A flexible distribution family for testing MCMC implementations

Papp, Tamás K. Working Paper A flexible distribution family for testing MCMC implementations IHS Working Paper, No. 60 Provided in Cooperation with: Institute for Advanced Studies (IHS), Vienna Suggested Citation: Papp, Tamás K. (2025) : A flexible distribution family for testing MCMC implementations, IHS Working Paper, No. 60, Institut für Höhere Studien - Institute for Advanced Studies (IHS), Vienna, https://irihs.ihs.ac.at/id/eprint/7280/ This Version is available at: https://hdl.handle.net/10419/324755 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ IHS Working Paper 60 August 2025 A flexible distribution family for testing MCMC implementations Tamás K. Papp All IHS Working Papers are available online: https://irihs.ihs.ac.at/view/ihs_series/ser=5Fihswps.html This paper is available for download without charge at: https://irihs.ihs.ac.at/id/eprint/7280/ Author(s) Tamás K. Papp Editor(s) Robert M. Kunst Title A flexible distribution family for testing MCMC implementations Institut für Höhere Studien - Institute for Advanced Studies (IHS) Josefstädter Straße 39, A-1080 Wien T +43 1 59991-0 www.ihs.ac.at ZVR: 066207973 License This work is licensed under the Creative Commons: Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/) All contents are without guarantee. Any liability of the contributors of the IHS from the content of this work is excluded. A flexible distribution family for testing MCMC implementations1 Tamás K. Papp ([email protected]) Institute for Advanced Studies, Vienna August 2, 2025 Abstract We propose a flexible, extensible family of distributions for testing Markov Chain Monte Carlo implementations. Distributionsarecreated bynestingsimpletransformations,whichallowvarious shapes,includingmultiplemodesandfattails. Theresultingdistributionscanbesampledwithhigh precision using quasi-random sequences, and have closed form (log) density and gradient at each point, making it possible to test gradient-based samplers without automatic differentiation. 1 Introduction Bayesian inference results in a posterior distribution over the parameter space of the model. Except for special cases, such posterior distributions do not conform to any well-known family, and their analysis requires sampling, which is a computationally intensive task that requires specialized algorithms and software developed for this purpose. The history of Bayesian computation is a procession of increasingly complex and sophisticated algorithms (Robert and Casella 2011), with constant new developments (Bou-Rabee, Carpenter, and Marsden 2024; Bou-Rabee, Carpenter, Liu, et al. 2025). Testing software implementations of these algorithms is challenging for three reasons. First, their output is stochastic, which requires a statistical approach, as opposed to simply comparing that the implementation mapsknowninputs to known outputs (Whittaker 1997; Ševčíková etal. 2006). Second, as emphasized by Grosse and Duvenaud (2014), it is difficult to separate the problems with the mixing of the Markov chains from the flaws in the implementation. Third, the very quantities we want to compare are, in general cases, unknown: as Talts et al. (2020) put it, The most straightforward way to validate a computed posterior distribution is to compare computed expectations with the exact values. An immediate problem with this, however, is that we know the true posterior expectation values for only the simplest models. This paper addresses the last point, by providing a very general distribution family which allows 1. accurate sampling using low-discrepancy sequences, which can be compared to MCMC output, 2. closed form calculations for the gradient of the log posterior, allowing the use of gradient-based methods (Neal 2011; Hoffman and Gelman 2014) without automatic differentiation, making the software implementation easier. We review the related literature in Section 2. Section 3 introduces the model family, while Section 4 demonstrates one possible approach to unit testing MCMC software. Section 5 concludes. 2 Related literature Geweke (2004) relies on posterior distributions being (unscaled) products of a prior and a likelihood. Samples are obtained in two ways: 1. sample p(θ)from the prior, then sample p(x|θ)from the generative model, 2. given a (θ,x)from the previous step, update θusing an MCMC transition, then generate data from p(x|θ). The combination of these steps should preserve the joint distribution p(θ,x). 1This paper documents the internals and rationale of the Julia software library Papp and contributors (2025b). Tamás K.Pappacknowledgesthe support of the AustrianNationalBankJubileumsfondsProjekt 18847. I would liketothank Robert Kunst for comments. 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Figure 1: Sobol sequence in 2D, 1024 points. Both methods should yield samples from exactly the same distribution, and they can be compared using various statistics. Cook, Gelman, and Rubin (2006) suggest the following procedure. Use the prior p(θ)to generate a parameters θ0, then the generative model p(x|θ0)to generate a data sample x. Finally, obtain posterior samples θi, for i= 1,2,...,L, and calculate the quantile ˆq(θ0) = 1 L∑︁L i=11θ0>θi. The authors show that as L→∞,ˆq(θ0)converges to the uniform distribution. This can be evaluated graphically or using formal statistics. Talts et al. (2020) suggest comparing the data averaged posterior pD(θ)=∫︂p(θ|y)p(y|θ)d(y)p(θ) to the prior distribution, where the integral is evaluated numerically, direct sampling fromt the prior p(θ)andthe data generating model p(y|θ), and MCMC forp(θ|y). Thecomparisonisperformed using one-dimensional rank statistics f:Θ→R, using empirical quantiles similarly to the ˆqof Cook, Gelman, andRubin(2006). Thepreferredmethodofcomparisonisvisual,andtheproposedsimulationstrategy also allows for error bands. 3 The distribution family It is assumed that the user has access to an implementation of low-discrepancy sequences2to generate quasi-random numbers u(i),i = 1,...,N in [0,1]m.3Figure 1 shows samples from a commonly used low-discrepancy sequence, the Sobol sequence, which we also use in our implementation. Each distribution implements the hypercube transform h:[0,1]m→Rnand a density function p:Rn→ R+that have the property ∫︂Rn g(x)p(x)dx≈1 N N ∑︂ i=1 g(︁h(u(i)))︁)ℓ(︁h(u(i)))︁ where the properties of the approximation depend on the variation of gas described in Morokoff and Caflisch (1995) and Caflisch (1998). Provided that a distribution has this property, we are able to 2See Niederreiter (1988). 3Lacking that, we assume that uiare obtained from a suitable random number generator, but that is less than ideal for accuracy so we suggest that it is avoided. 2 1. obtain low-discrepancy samples using h, and 2. compare them to the the result of MCMC methods sampling from f. We describe the construction of a family of functions that have this property below. For numerical accuracy, we construct the log density function ℓ(x)=log(p(x)) and its gradient ∇ℓ(x). Notice the two dimensions associated with a distribution in this family: nis the dimension of the domain, while m≥nis the hypercube dimension. The two may be different because we need extra values to generate mixtures, as explained in Section 3.3. 3.1 Primitives Assume that a univariate distribtion Dis available with log density ℓ:R→R+and inverse cumulative distribution function I:[0,1]→R. Then, for a given dimension n≥1, we define ℓD(x)= N ∑︂ i=1 ℓ(xi) ∇ℓD(x)= N ∑︂ i=1 ℓ(xi) hD(u)=⎛ ⎝ I(u1) ... I(un)⎞ ⎠ Averyconvenientprimitivefamilyisthestandardnormal. WeintroduceStdNormal(n)usingtheabove notation4with ℓ(x)=−x2/2+log(2π) 2 I(u)=√2erf−1(2u−1) It follows that ℓStdNormal(x)=−∥x∥2 2+nlog(2π) 2 ∇ℓStdNormal(x)=x hStdNormal(u)=√2⎛ ⎝ erf−1(2u1−1) ... erf−1(2un−1)⎞ ⎠ Figure 2 shows the contour plot of the highest density regions of the StdNormal(2) distribution, primarily for establishing the convention of illustrating distributions in R2with 90%, …, 10% highest density contours in the paper.5 3.2 Transformations All transformations are defined using a differentiable bijection g:Rn→Rnwhere nis the dimension of the distribution. In this section we use the convention y=g(x). 4Notation in this paper closely follows the software package Papp and contributors (2025b), with names exported from that package typeset in slanted font. Upper/lowercase and minor spelling conventions may differ. 5Given the exact sampling using Sobol sequences, such contours are trivial to obtain computationally. 3 −3−2−1 0 1 2 3 −3 −2 −1 0 1 2 3 Figure 2: The standard normal distribution StdNormal(2). Contour lines are highest density regions containing 90%, …, 10% of the probability mass. Given a source distribution S, characterized by ℓS,hS, the destination distribution Dhas hD(u)=g(hS(u)) ℓD(y)=ℓS(︁g−1(y))︁−log(︃det(︃∂ ∂x(︁g−1(y))︁)︃)︃ ⏞ ⏟⏟ ⏞ ≡c ∇ℓD(y)=∇ℓS(︁g−1(y))︁∂ ∂y (︁g−1(y))︁ ⏞ ⏟⏟ ⏞ ∂x ∂y −∂c ∂y where cis the correction required for transformation of the variables. For all transformations, we provide analyitical characterizations of g−1,c,∂x ∂y , and ∂c ∂y , which makes calculations possible without automatic differentiation. All transformations preserve the hypercube dimension of the source distribution. 3.2.1 Linear map The linear map Linear(A)uses the mapping y=g(x)=Ax, for an n×ninvertible matrix A. Then x=A−1y∂x ∂y =A−1c=log(|det(A)|)∂c ∂y =0 Figure 3 shows an example. 3.2.2 Translation Translation Shift(b)takes a vector bof length n, and is defined using y=g(x)=x+b. Then x=y−b∂x ∂y =0 c=0 ∂c ∂y =0 Shift may seem like a trivial transformation, but it is very useful since it allows omitting a location parameter from the other transformations: a general transformation τcentered around the origin can be recentered around an arbitrary location bas Shift(b)◦τShift(−b). We use this extensively below. 4 −3−2−1 0 1 2 3 −3 −2 −1 0 1 2 3 Figure 3: Linear(A)(StdNormal(2)) with A=(︁0.2 0.5 0.4−0.7)︁. Contour lines are highest density regions containing 90%, …, 10% of the probability mass. 3.2.3 Elongate TheElongate(k)transformationstretches(orshrinks)thetailsaroundtheorigin. Foragivenparameter k∈R, it is defined by y=g(x)=x·(︂1+∥x∥2 2)︂k In order to characterize g−1, we define ξ(Y,k)as the X >0that Y=X·(1+X2)k(1) for all Y > 0. The solution always exists since the right hand side of (1) is increasing in X. Calculating ξ(Y,k)requires a univariate numerical solver, such as bisection or Newton’s method.6 Then, for a given y, let κ=ξ(∥y∥2,k)2D=(1+κ)−kx=D·y c=knlog(1+κ)2+log(1+kκ/(1+κ)) A=1+κ B =1+(1+2k)κ Using these, we can calculate ∂x ∂y =(In−(2k/B)xx′)D∂c ∂y =k(2+Bn) A2k+B2·y Figure 4 shows an example. 3.2.4 Funnel The transformation for Funnel, inspired by the well-known example of Neal (2003), uses the mapping y=g(x)=⎛ ⎜ ⎜ ⎝ x1 x2exp(x1) ... xnexp(x1) ⎞ ⎟ ⎟ ⎠ 6Newton’s method can be set up in a way that it is always convergent in a few steps. See the related source code. 5 −3−2−1 0 1 2 3 −3 −2 −1 0 1 2 3 Figure 4: Elongate(0.5)(StdNormal(2)). Contour lines are highest density regions containing 90%, …, 10% of the probability mass. −3−2−1 0 1 2 3 −3 −2 −1 0 1 2 3 Figure 5: Funnel(StdNormal(2)). Contour lines are highest density regions containing 90%, …, 10% of the probability mass. which has ∂x ∂y i,j =⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1if i=j=1 exp(−y1)if i=j=1 −yiexp(−y1)if i=1,j =1 0otherwise c=(n−1)y1 ∂c ∂y =⎛ ⎜ ⎜ ⎝ n−1 0 ... 0 ⎞ ⎟ ⎟ ⎠ Figure 5 shows an example. 6