I started thinking a bit more about what we’re doing when we use something like posteriordb (from Stan) or Inference Gym (from TensorFlow Probability) to evaluate a sampler.
posteriordb gives you Stan programs and 10,000 reference draws from their posterior. The reference draws must be generated by another sampler. In some cases, like Neal’s funnel or high-dimensional normals, etc., we can just run forward simulation and get independent draws. In cases of models with data, we don’t actually know the posterior analytically—that’s why we’re running MCMC in the first place! I think what posteriordb did was run NUTS for enough iterations it could be thinned down to 10,000 roughly independent draws.
Standard error in reference
The standard error for $latex N$ independent draws estimating a particular quantity like the posterior mean of a parameter is
$latex \textrm{se} = \frac{\displaystyle \textrm{sd}}{\displaystyle \sqrt{N}}$
where sd is the standard deviation of the quantity in the posterior. Standard Monte Carlo methods scale horribly due to that square root—Markov chain Monte Carlo can be even worse because rather than the number of draws, we have the effective sample size, which can be much lower than the number of iterations for more challenging posteriors due to the correlation in the draws.
Because $latex \sqrt{10,000} = 100$, with 10,000 draws, or standard error is only 1/100th of our standard deviation. In general, we need four times as many draws to cut the error in half, or 100 times as many draws to add a decimal place of accuracy by cutting error by a factor of ten. That means we can’t reasonably expect something like posteriordb to store enough draws to make fine-grained convergence evaluation possible. A slightly more economical alternative would be to store long term averages. There we cold run a million iterations and save parameter and squared parameter estimates. We could even run a hundred million by farming out 100 jobs on a cluster. But that’s still only three or four decimal places of accuracy. After that, maybe we could wring out another decimal place or two with control variates.
That’s a lot of residual error to have to deal with when evaluating an algorithm that asymptotically approaches the correct answer. The only thing that saves us is that the algorithms we’re evaluating are also samplers and are subject to the same inverse square root rate of reduction in standard error.
Report error in evaluations
What we really want to report is the standard error of the residual between our new model’s estimate and the truth. But we don’t usually have the truth. In practice, we usually do one of two things, depending on how compute-intensive it is to fit my model and our compute resources and patience.
Easy model:. Run the new algorithm 100 times or even 1000 times and report something like a histogram of errors. With this many runs, we can directly estimate standard error (assuming our draws are unbiased) and we can estimate it even more tightly and without assumptions if we know the true values.
Hard model:. Run the new algorithm once or twice and use an estimated standard error for the new algorithm and report that instead of an empirical histogram of errors. This suffers from estimation error in the effective sample size, which can be large.
The problem with both of these aproaches is that it’s only reporting variability of the new algorithm, treating the reference draws as true. We should be going one step further and folding the error in posteriordb or another reference database into the uncertainty reporting of the error of the new algorithm. The additional error could be done with simulation (and maybe analytically—normals are nice that way).
That’s easy enough to do reporting single univariate means. I do not know how to do this for evaluating error in something like KL-divergence, which itself has a Monte Carlo integration, or something like Wasserstein distance. In the Pathfinder paper, we used a semi-discrete Wassertein distance algorithm. I have to admit I don’t even know how that algorithm works.
Evaluating delayed rejection
I ran into this problem because I want to evaluate our delayed rejection HMC and show that the error goes to zero as I take more draws. Specifically, I want to extend the illustration in the paper that shows that standard HMC and NUTS fails, but DR-HMC succeeds at sampling from multiscale distributions like Neal’s funnel in its centered parameterization.
There are two fundamental problems with this.
Problem One: By the time we get past an effective samle size of 10,000 from the new algorithm, the error we’re estimating is dominated by the standard error of posteriordb. We don’t even know the real error of posteriordb, only an estimated standard error assuming draws are independent.
Problem Two: The second problem is that we have no guarantees that the reference draws are unbiased draws from the posterior. For example, in the centered parameterization of Neal’s funnel, there is no step size for HMC that will properly explore the posterior—each step size picks out a range of log scale parameters it can explore (as we demonstrate in the delayed rejection paper). The standard error estimate from the reference draws assumes independent draws from the posterior, not biased draws. If the posteriordb draws are biased, the error reported for a perfect sampler will not go to zero as the number of posteriordb and new sampler draws goes to infinity.
Help?
Does anyone have any suggestions of what we can do in the face of potentially biased and stochastic samples for evaluation? I could wait until the Stan meeting on Thursday and ask Aki, who’s a font of inspiration on topics like this, but I thought I’d ask a wider group and see how that goes.
Bob, as you were an author of;
Delayed rejection Hamiltonian Monte Carlo for sampling multiscale distributions
Chirag Modi
Alex Barnett
Bob Carpenter
https://www.researchgate.net/publication/355061236_Delayed_rejection_Hamiltonian_Monte_Carlo_for_sampling_multiscale_distributions
I’d suggest asking Terry Tao.
This reminds me of an experiment I haven’t got around to doing, but should be pretty simple. Start with a newtonian solar system simulation using timesteps of say 1 day. Obviously this will deviate from observations.
But, the idea is now take that simulation as the benchmark. Then run the same simulation using a timesteps from one year to one nanosecond (or whatever is feasible). Won’t the “error” (with respect to the chosen 1 day benchmark) be minimal at 1 day and grow larger as you get more or less precise from there?
Same with using analytic solutions as the benchmark. If there really is no infinitely small unit of time your totally correct model will start getting deviation from observation in cases where those very small (~10^-44 sec or whatever) errors can accumulate. Ie, over vast distances or time. The problem would be with the implementation, not the model itself.
In your case I’d say you should be optimizing for not only goodness of fit (to the, as you say, imprecise benchmarks) but also how much it costs to get there (how long it takes, and what resources are required). So the “optimal” sampler will depend on the specific circumstances/requirements.
tldr: It doesn’t make sense to have a single best way of doing things. It depends on the problem you want to solve along with time/resource constraints.
Not really. That’s precisely why stan uses a symplectic integration scheme
Perhaps I need to be more explicit. I was talking about comparing solar system simulations to a benchmark, *in analogy* to comparing whatever sampler in stan to a benchmark.
Whether stan is symplectic is not relevant at all.
And anyway stepsize is a key parameter for symplectic integration, which is why it is usually varied to control the error. Pretty sure if you ran that experiment the fit would be best when using the same timestep (in fact it should be exact if the model is deterministic), even using symplectic integrators. But, I also want to actually do it to see the curve.
Yes, that’s right. What the leapfrog integrator does is sample a “shadow Hamiltonian” exactly, so you can get a grip theoretically that way. The shadow Hamiltonian varies with step size and approaches the true Hamiltonian as step size goes to zero and leaves stability beyond a size determined by the spectrum of the Hessian. Having said that, leapfrog isn’t a very accurate ODE solver per se. What it is great at doing is preserving medium scale error and preserving a quadratic error long term. The long-term stability is great for Metropolis proposals, though any work beyond half an orbit or so is wasted (though often required by detailed balance in existing algorithms, like where NUTS integrates in both directions in time).
Varying step size in such a way as to adapt to curvature and preserve the stationary distribution in MCMC is tricky. We can do it with delayed rejection and we’re working on a second alternative that’s similar in spirit if not in details to the recent autoMALA paper.
I use the following to test my own Julia implementation of NUTS (DynamicHMC.jl):
1. primitive distributions which you can sample from exactly, using low-discrepancy sequences (basically, anything composed of independent univariate distributions). I use Sobol sequences.
2. a bunch of transformations, including affine, a nonlinear heavy-tail one (y = (1 + ‖x‖²)ᵏ⋅x), and mixtures. They propagate the mechanics to sample using Sobol sequences.
Nesting these transformations, you can compose challenging distributions for testing, and preserve the property that you can sample from them exactly. Maybe I could add a funnel-like one too.
The code is in https://github.com/tpapp/LogDensityTestSuite.jl, PM me if you have questions.
Tamas:
Also, in simple hierarchical models you can integrate out the linear parameters in the posterior distribution, leading to a low-dimensional posterior that can be computed to arbitrary accuracy using a grid search, as illustrated for two different examples in chapter 5 of BDA. Bob must be aware of this regarding the funnel (or, as he puts it, “Neal’s funnel”), so I guess he’s talking about the challenge of checking more complicated problems.
Tamas—I’m not that familiar Sobol sequences or with quasi-Monte Carlo, which is the only place I’ve seen them used. I’m not clear on what you mean by sampling from a distribution exactly. Does that mean you can take a sample of N i.i.d. random draws? I thought the point of Sobol sequences is that they are deterministic and designed to be spread out, not independent, so that expectations converge faster than 1 / sqrt(N). They’re certainly atypical in the way they’re spread out. Also, aren’t they deterministic? What prevents the next sample of N draws from being identical?
Andrew—I thought everyone called it “Neal’s funnel” because he introduced the density. Yes, we can sample exactly from Neal’s funnel, because it’s a simple forward model, which is where we can get reference draws. We know the means are all zero and know the marginal variance of the log scale parameter, but I don’t know how to analytically evaluate the variance of the coefficient parameters. Hence the blog post.
When the models get more complicated and data gets involved, we don’t know how to integrate those parameters out in general. We can do it in some simple normal-normal cases, or simple and low enough dimensional cases they fit in BDA, but we can’t do it in general or Andrew wouldn’t still be asking the Stan team to make hierarchical models go faster :-). INLA does this kind of marginalization with painstaking attention to detail on a model-by-model basis. Hopefully Charles Margossian’s black box method will work out—he’s too busy doing research to get a job to have much time to work on it, though—academic hiring committees tend not to care if you solve practical problems for tens of thousands of people.
I think the point of the Sobol sequence is that it’s both fast to compute and does converge rapidly, so you can estimate quantities with precision efficiently. So you could draw say 1000 samples via NUTS and compare it to 10000 “exact” Sobol samples (which can be drawn in the blink of an eye) and the Sobol samples statistics could be considered numerically close enough to exact that you could detect meaningful discrepancies in the NUTS samples. Ie. its a numerical technique for getting something fairly numerically precise.
Bob:
The first fully Bayesian analysis of a normal-normal hierarchical model is in Tiao and Tan (1965). Or, at least, that’s the first one I know about. Maybe Laplace or Gauss did it too somewhere.
Radford Neal is wonderful, so though he’s not the first to work with the funnel, I have no problem with him getting some credit in the world.
Bob, Andrew: Just to clarify, I am not proposing to integrate out anything.
What I am using is a family of distributions all transformed from a [0,1]^N hypercube, which is very easy to sample from using Sobol (or any other low-discrepancy) sequences.
The idea is that nesting simple nonlinear transformations (incl mixing) can yield quite challenging distributions (for NUTS and in general), which have a (log) pdf that is differentiable and at the same time an exact sampling method because they eventually come from a hypercube.
The challenge is trying to mimic distributional properties that stress-test MCMC algorithms within this family. So eg when I debug a sampling problem and identify a challenge, I need to boil down the essence of it to a transformation.
I’m seconding this approach. At their best, low-discrepancy sequences converge proportional to 1/N, a huge improvement over 1/sqrt(N). The tricky part is dealing with multi-parameter sampling, as it is easy to contaminate your random draws with cross-parameter correlations that come from the sequences themselves or reintroduce randomization and revert to 1/sqrt(N) convergence. This bit of C++ code seems relevant (the code itself is on the very bottom of the page):
https://people.sc.fsu.edu/~jburkardt/cpp_src/sobol/sobol.html