Nawaf Bou-Rabee was telling me and Chirag Modi that a high-dimensional normal initialized at the mode would be a disaster. I was skeptical as I’d never seen this behavior. So I dashed back to the computer and ran some simulations. Summary: Nawaf’s right! And the reason I hadn’t seen it is because NUTS handles the problem with variable adaptation during warmup.
Edit: I forgot the firsrt time to say why after I had written it into the PyMC Discourse response, but Daniel Lakeland reminded me. It’s because the mode is highest density, so all your proposals will be lower density. And the density goes down quickly in high dimensions as you move away from the mode. I have some illustrations of this in my case study on the curse of dimensionality.
Simulation results
I defined a 10,000-dimensional standard normal in Stan and used CmdStanPy to run Stan’s default NUTS for 4 chains of 1000 iterations each to give me 4000 draws.
If I initialize at the origin and fix the step size to 1 (which is stable for the leapfrog algorithm), all of the transitions diverge. I get an effective sample size (ESS) of 0.
If I initialize at the origin and fix the step size to 0.5, I get an ESS of about 20.
If I initialize randomly uniform(-2, 2) using the Stan default and fix the step size to 0.5, I get an ESS of about 2000 (in 4000 draws).
If I initialze at the origin, but let NUTS do its thing adapting step size, I get an ESS of about 9000. NUTS can use multiple step sizes during adaptation, so it can get away from the origin with a small step size and then crank it back up during the last phase of adaptation before sampling.
If I take Stan’s defaults, which is to let NUTS do its own thing adapting step size and initializing randomly uniform(-2, 2), I also get an ESS of about 9000.
Cool! The defaults for NUTS are very robust even in this nasty example.
PyMC provides this option and used to recommend it, but they must have done this ages ago and they now discourage initialization at the mode in their documentation.
Stan program
parameters {
vector[10000] y;
}
model {
y ~ normal(0, 1);
}
Python script
import cmdstanpy as csp
import numpy as np
N = 10000
y = np.zeros(N)
model = csp.CmdStanModel(stan_file='normal.stan')
# INIT ORIGIN, STEP SIZE >> 0.5
fit0 = model.sample(inits={'y': y}, step_size=1, iter_warmup=0, adapt_engaged=False)
fit0.summary()
# INIT ORIGIN, STEP SIZE = 0.5
fit1 = model.sample(inits={'y': y}, step_size=0.5, iter_warmup=0, adapt_engaged=False)
fit1.summary()
# INIT RANDOM, STEP SIZE = 0.5
fit2 = model.sample(step_size=0.5, iter_warmup=0, adapt_engaged=False)
fit2.summary()
# INIT ORIGIN, ADAPT STEP SIZE
fit3 = model.sample(inits={'y': y})
fit3.summary()
# NUTS DEFAULT INIT (uniform(-2, 2)), ADAPT STEP SIZE
fit4 = model.sample()
fit4.summary()
Intuitively, once you’re at the mode, every step you try causes the log posterior density to decline drastically. So you’re going to reject the transition with high probability.
On the other hand with HMC initialized away from the mode you are roughly in the position of a skateboarder at a skate park up high on the walls around the edges… when you step into the skate park you will fly “down” the hill toward the mode, so you’ll move in the right direction. The main thing is you’ll potentially get to very high velocities if you’re “very high up the wall” and this will become unstable without small timesteps.
Once your chain has converged, in a high dimensional space, each new sample will be at essentially precisely the same log density (plus or minus a small fluctuation). The *ideal* place to initialize your chain is at any point which is roughly within the [0.01,0.99] quantiles of the log density. If you do that, you are essentially already in equilibrium or close to it. The trick is to estimate the range of the log_pdf while finding the initialization point.
Thanks, Daniel. I forgot to mention why it’s a problem! I had just responded to the PyMC list and put the reasoning there, but forgot to include it here.
Initializing from the stationary distribution is ideal as then you’re stationary to start and don’t need any burnin to reach approximate stationarity. You do that by taking a draw with a log density in a central interval of posterior log densities. The typical set results in information theory are about “typical” draws having log density near the expected value. That was the whole motivation for the Pathfinder algorithm—run an optimization algorithm and figure out when the log density is close to its expected value. We tried a few ways to do that (runnig short chain MCMC and estimating volume times density) before turning it around and thinking about it like a variational inference problem.
Bob:
Just a minor correction here regarding your use of the term “typical set.” You are using this term in the way that we have often used it when discussing HMC. But this is not quite the same as the meaning of the term in information theory.
From my earlier post on the topic (which was heavily informed by conversations with you):
Here’s my summary of the differences between the official definition of “typical set” and some of how it’s discussed in the Stan community:
As I wrote in comments, I think the Stan community might be better off abandoning the term “typical set,” at least we can agree on a definition. Talking about log p rather than the typical set has been very liberating.
There’s lots more at the above-linked post and its comment thread.
Andrew I don’t disagree with anything you say, but I just wanted to mention that the information theory motivation is basically observing many repeated samples of some 1 dimensional measurement (like, say heights of people). But if you are sampling from a model which is an independent 10 dimensional normal then 1 sample from that 10 dimensional vector could also be thought of as 10 samples from 1 dimensional normals (so each single sample *is* a sequence of N).
Where it gets tricky is that in Bayes there’s some dependency between dimensions so it’s not quite like the sequence of independent measurements. So in some sense our single sample of an N dimensional models in bayesian stats behaves more like an N dimensional sequence of data values than it does a single 1 dimensional measure (like 183 cm) but in other respects (dependency) it doesn’t quite. The concepts that apply to sequences of N measurements also apply pretty well to single samples of parameters in N dimensional models.
If you said “the typical set is the set of points x where log(p(x)) is greater than p*” where p* is appropriately defined then the sampler I mention above, the “white box sampler” is a sampler that samples the distribution which is proportional to the indicator function on that set.
Once you’ve sampled the indicator on that set, getting it distorted by the actual density is relatively easy compared to exploring the set in the first place. Hence the motivation.
Daniel:
I guess the main point of my post was that the casual use of the term “typical set” was leading to confusion, in part because there is this thing called “typical set” that has a technical meaning, but that’s not what the Stan people were talking about when they used that term.
Often this seems to happen in statistics, that an expression will be used that has a sort of misplaced precision, and then misunderstandings arise.
I really did mean to write “typical set” in the technical information-theoretic sense where an epsilon typical set defined by
T_e = { theta : log p(theta) in E[log p(Theta)] +/= epsilon }.
I find the actual definnition Daniel used of central 98% of log densities easier to understand for stats purposes,
C_alpha = { theta : theta in central alpha interval of log p(Theta) }
so it’s what we used to measure convergence in Pathfinder: how quickly could we get into C_0.99 from a random init.
If you know the mode and have an estimate of the curvature (Hessian), does it make any sense to initialize based on random draws from the implicit, approximate multivariate Normal posterior that you would assume based on those quantities? (I know, I could just try it for myself and see …)
Basically draws from a variational inference posterior, yeah, it does make sense to initialize that way. Pathfinder is probably even better.
I’ve got an algorithm I would really like to implement for use in Julia, but I haven’t ever gotten clear instructions on how to write a sampler for Turing. I think I could probably figure out one for generic LogDensityProblems though. The methodology is to start near the mode, and subtract some fixed quantity from the modal log-density and run a sampler that goes forward in a random straight line and simply reflects off the boundary (toward the interior in a random direction) when it crosses the lp value mode – fixed_quantity. You can decide if it’s towards the interior by taking a step and seeing if the lp increases or decreases if it increases continue if it decreases change directions spherically at random until it increases.
Sample that process at uniform intervals in time (it’s a piecewise deterministic markov chain) and collect that sample. Then after that sample is collected, take each point in the sample in turn (Gibbs style) and use diffusive Metropolis-Hastings MCMC for some fixed number of steps on each point until the sum of the log densities becomes constant, and call it a day.
There are basically three tuning parameters here, one is the initial offset of the log-probability that forms the boundary, and one is the timestep size for the piecewise deterministic model, and one is the scale for the MH-MCMC (basically use the initial sample to find correlation structure and then downscale it for the MH steps.
Note that nothing about this sampler requires derivatives at all. The reflection process avoid derivative calculations by simply requiring regeneration. It could get expensive if you’re stuck in a weird funnel or whatever, but for many spaces it’s gonna work well (I’ve tried it). I call it a “white box sampler” because the piecewise deterministic process is basically the same process as a photon inside a box painted with 100% reflective but white (not silver/mirror) paint.
You can think of the white-box sampler as an “umbrella” sampler that samples from the “uniform within the boundary” space. This space is essentially exactly the same as the actual space for higher dimensional models (ie. 10+ certainly by the time you’re in the 50-100+ range). The final metropolization step is just there to distort the “almost perfect” sample into a “perfect” sample.
The big benefit is that no derivatives are required.
Edward Roualdes’s motivation for developing BridgeStan initially was so that he could develop algorithms in Julia using Stan models. We use it for everything algorithm development related now. We have few enough users there that we’re happy to help with problems that come up.
The algorithm you’re describing rings more than one bell.
The Apogee-to-Apogee Path Sampler (Sherlock, Urbas, Ludkin 2022) is an alternative to NUTS that uses U-turns in target log density (which is one-dimensional) rather than in parameter space. The authors found more than one U-turn was better for correlated posteriors with a lot of transverse movement around the principal direction and then needed to do the same kind of trick as NUTS forward and backward in time to maintain detailed balance. It uses gradients, of course.
The searches you describe can get tricky numerically. One way to do it is to use slice sampling over the paths picked out by the Goodman and Weare ensemble sampler. That adjusts for global scale and curvature with a few global draws.
The quantum physics folks really really like to lay down a uniform hypercube and then try to transform to the target density they want to sample. It’s somewhat related to quasi Monte Carlo. You can even think of something like normalizing flows the same way—transform a unit hypercube to a good variational approximation of the target density; a Laplace approximation where available can be a better base distribution than a hypercube for faster and more stable convergence.
Your idea also reminds me of Repelling/Attracting MCMC (Vishwanath, Tak). It alternates negative and positive friction to move through log densities faster. There are some cool multimodal applications in the paper.
So yes, I think there’s a lot to explore here.
You may remember that we also tried this sampler out on the mc-stan discourse 4 years back:
https://discourse.mc-stan.org/t/the-typical-set-and-its-relevance-to-bayesian-computation/17174/19?page=2
it may have been one of the early nudges towards working something that eventually led to Pathfinder, not sure.
The Julia implementation I had back then was a skeleton. if I can figure out how to write a sampler for a LogDensityProblems.jl object then maybe I can make it more general and we can test it against other models in the posteriordb or something.
Ben:
That sort of mode-based normal approximation can work for simple problems but not for hierarchical models, which is one reason we developed Pathfinder.
@Ben: There are two problems with using the Hessian in a Laplace approximation. (1) Hierarchical models don’t have a mode, and (2) quadratic is really expensive in high dimensions. That’s why we developed Pathfinder, which tries to find a posterior mean and uses a low-rank plus diagonal approximation to keep costs linear.
Thanks for this post! Anecdotally, I’ve seen this behavior for a long time but never understood why. Your and Nawaf’s explanation (and Daniel’s great metaphor) makes it very clear. Speaking for myself, this makes me realize that I’d imported practices and intuitions from my earlier experience with numerical optimization (and from BUGS/JAGS) that simply don’t apply to HMC.
That’s why I wanted to post. I hadn’t thought about this the right way until Nawaf pointed it out. It’s also a problem for random-walk Metropolis. Either way, you’re trying to climb an energy barrier if you think of the log density as a potential.
*negative log density as a potential (high energy is “hard to get to” is low probability)
This lead me to the following paper:
https://pubmed.ncbi.nlm.nih.gov/11384462/
If you manage to initialize at the mode you already have the exactly correct answer right? So ideally what you want is for all proposals to be rejected right away. This means you saved a bunch of computations.
Maybe it can be seen as a model misspecification issue. You don’t believe there is an exactly correct answer due to whatever uncertainty, but that is missing from your model. Perhaps related to the equilibrium-based physics motivation, which is more like an average over a discrete set of instantaneous states.
That appears to be a strategy for finding a “ground state” which is to say the optimal lowest energy point. so if you start there, yeah you’re done.
But it doesn’t seem like it’s a strategy for general sampling from a distribution, just some kind of optimization method?
Yea, its essentially simulated annealing except there is a discrete threshold with no random Metropolis-Hastings step. The adaptive phase of MCMC amounts to using optimization to initialize the values. But look what happens when you already optimized “too much”:
It blows up. Isn’t this because you started with the right answer which is essentially infinitely better than anything else that could ever be proposed?
infinitely better? no. Notice that he can make it work by dropping the step-size so it can crawl out of the hole.
>If I initialize at the origin and fix the step size to 0.5, I get an ESS of about 20.
What happens if you use stepsize of 0.5 and sample from:
y ~ normal(0, 0.5);
Right, if you make the step size infinitesimal then it will be stable for any distribution, as we increase the step size to a standard positive value eventually there’s a stepsize that is unstable, and the size of that depends on the distribution you’re sampling from.
Shouldn’t that be a function of whatever process generated the data along with the model you are fitting to it? It wouldn’t normally show up since typically all models are misspecified in one way or another. But if everything exactly matches up I bet that is when you see this happen.
That shouldn’t be a bad thing though. How is it a disaster to not waste a bunch of time and energy when you already have the right answer?
Anon, sort of, the step size instability is a numerical computing issue with a sampler, that is, a process that moves around the parameter space and spends time in regions of high probability. Similar to solving a differential equation can become unstable when you take a time step too long, sampling from a posterior distribution can become unstable when you take a time-step too long in an HMC sampler (because an HMC sampler is in fact a solver for a differential equation).
The only thing it has to do with the model and data is that the shape of the curves through space that the HMC sampler follows can be quite convoluted depending on the model and the data. The more curvature of the paths the smaller the stepsize will need to be to remain stable. It’s not really a modeling issue, it’s a sampling/numerical issue (sort of, models that are simpler will tend to suffer less from this kind of numerical issue)
Say Bob changes his model to:
1) Generate data from normal(0, 0.5)
2) Use y~ normal(0, 0.5) as the model,
3) Uses stepsize of 0.5 (which I gather also amounts to the standard deviation of a normal distribution)
Will he get results just like he saw from using 1 as the standard deviation? Ie, it blows up and effective sample size ~0?
I didn’t try it.
Yes, unless I’m mistaken, this is completely symmetric with the problem of using normal(0,1) and 1 step size. When it comes to the convergence ability of a given step size the dimensionless ratio of step-size to standard deviation is all that matters.
Think about it this way, when you have a particular pdf curve, like a normal bell curve, how can you measure “how far out on the curve you are” ? The answer is you have to compare how far you are to some measure of the size of the curve. One such measure is the standard deviation, another such measure might be say the interquartile range. Another might be say the 90th percentile… You’ll find that all these measures are proportional to each other because the curve is defined in terms of the standard deviation. So when you’re measuring “how far do you go in 1 unit of time” you do it in terms of “fractions of the standard deviation”. So normal(0,0.5) and 0.5 step size is the same as normal(0,1) and 1.0 step size since the step size to standard deviation ratios are equal.
Makes sense. I still didn’t play around in stan to make sure…
So if you:
1) initialize/find the maximum density
2) data-generating distribution == modeled distribution == proposal distribution.
3) Something like max(density)/expected(density) > ???
Then you get ESS = 0. Ie, not just low ESS. I think I failed to make clear that this is the case I am concerned with.
And the opinion of everyone is that this situation cannot be handled gracefully, even in principle. You need hacks like “don’t initialize to the mode” or the elaborate ad hoc step-size adaptation we all use.
ESS=0 comes about essentially because every time you try to take a step using the large step size and starting from the mode you reject the step, so you just get a markov chain that’s the same point over and over and over again.
With a small step size you can walk around and eventually make your way into the right region but it’s inefficient. if you started in the right region you’d be able to take somewhat larger step sizes the whole time and more rapidly get a useful sample.
No.
1. You can make any point you want into the mode with an arbitrary coordinate transform. That doesn’t mean you can pick any point you want to be the right answer
2. Even if you can prove that you have a special parameterization where the expected value of the parameters is near the mode, that high concentration doesn’t mean that the variation doesn’t matter for practical purposes. For a 10,000 dimensional unit normal, the expected absolute deviation from the origin is still 10,000 * sqrt(2/pi) even though the density is extremely high at the origin and falls off rapidly. To make this practical, any betting on the individual dimensions still has to take into account the full unit variance no matter how sharply peaked the joint density is. That’s not “missing from the model”.
I guess it depends on the question. The paper Anon links to is about mode-finding I thin. If your question is “where is the mode” then yes if you start at the mode you’ve started at the right answer to that question.
But for a sampler which aims to draw a set of values which together qualify as effectively as good information wise as some N of IID draws (ie. ESS = N) then you have a different problem and everything you said applies yes.
You can get the paper for free here:
https://www.researchgate.net/publication/28353178_Best_Possible_Strategy_for_Finding_Ground_States
It is about optimization, and not mapping the posterior. The relevant point was that it turns out using a discrete threshold works better than the Metropolis-Hastings approach of rarely accepting values with far worse fit. As quoted above, detailed balance does not hold in this case.
The problem in the OP is that it started with the best fit, and then the output of the MH approach has been interpreted as some kind of problem. Seems to me actually finding the right answer is discouraged by the algorithm used, due to the MH step.
I mean yeah, if your goal is optimization, then the mode is the right answer. The OP is obviously not about optimization though, so I’m not sure why you would bring it up. If you wanted to optimize, gradient descent and its descendants would generally be better than random walk methods anyways; the reason to use random walk methods at all is if differentiation is not straightforward, in which case a tool like stan or HMC isn’t going to work for you anyways.
It really just seems like you don’t understand the difference between optimization and integration. The solution to the 10,000 dimensional normal optimization problem is the zero vector. How does knowing that single value help you if you were, say, betting on the length of one draw? If optimization were always the “correct answer”, there wouldn’t be any need for monte carlo algorithms at all
Afaict bob asked the computer “What is the mean of these 10,000 distributions I simulated from normal(0, 1)?”
In that case, the correct answer is exactly zero for all of them.
He didn’t ask for the predictive distribution (ie, what would happen if you took samples). The standard deviation of the posterior corresponds to SEM, not sd.
The computer said you already gave me the answer so SEM = 0, and I’ll output something that looks like nonsense to you because you expected some uncertainty.
1. The mean is not generally the mode. If this were a skewed distribution, the mode would not be the mean but initializing at the mode would fail in the same way; optimization would not give the “correct answer” to the question of the mean
2. We’re not necessarily asking fit the mean anyways; the samples can let you estimate the expectation value of an arbitrary function of the draws. Eg you could average the norm of the draws to estimate expected distance from the origin, which no single value will help you with
That’s not what he asked for. I thought I wrote something here, but apparently I didn’t. But briefly, how would you compute the expected distance from the origin with your one ”correct answer”? And what if the mode were not the mean
Anon, what Bob wanted was a monte carlo sample from the distribution. If he starts from the mode he has difficulty moving the point around, but the point of a monte carlo sample is to give you a sample such that 1/N*sum(f(x_i)) is close to the mean of f(x) for every function f
This must be Bayes bait! Most of our models don’t have modes—they’re hierarchical (or what a frequentist would call “mixed effect”). Plus, we don’t just want the mode, we want to be able to estimate expectations. We do this for parameter estimates as posterior means (that minimizes expected square error with finite data whereas modes only have asymptotic convergence guarantees), for predictive inference and for event probability estimation.
Thanks for the response. I understand the intent of MCMC. The problem is that the code is not doing what you intended because (from your pymc link above):
This essentially means you managed to find the optimal model (including whatever parameters). Then what is the optimal proposal distribution? I found:
https://ucla-biostats-285.github.io/reading/AdaptiveMCMC.pdf
So when everything is perfect and optimal, the algorithm doesn’t behave as intended. Its like the dog that finally catches a squirrel then just looks around confused.
Oh dear god every second anyone has spent trying to communicate with you has been a waste of time
Anon,
The goal of MCMC is to take samples from the region of space where the log(p(x)) is “close to” integrate(log(p(x)) p(x) dx) (that is, close to the expected value)… Because almost all the points you’d get from a “perfect” sampler would be in that region.
usually at the mode log(p(x)) is much different from the expected value, so you *want* to get away from the mode, but the MCMC mechanism has a hard time getting away from that region because the log(p) declines rapidly and the difference in log(p) tells you whether you should accept the transition. In essence you’re stuck in a region of space that is very unusual.
By setting the timestep very low you can walk your way away from the mode into the proper region of space eventually, but it’s not efficient, it takes a lot of computation. You’d be much much better off starting away from the mode and closer to the region with the desired log(p) because you’ll get into the desired region of space more quickly and use larger timesteps and have lower correlation between consecutive samples.
Mode finding and MCMC are very different beasts.
Actually the “goal” of MCMC is not to take samples from the region of space where the log(p(x)) is “close to” integrate(log(p(x)) p(x) dx) any more than it is to take samples from the region of space where the log(p(x)) is “not lower than” than integrate(log(p(x)) p(x) dx). Because almost all the points you’d get from a “perfect” sampler would be in that region as well.
Carlos, if you say “not lower than mean(log(p))-epsilon” for a sufficiently large epsilon then yes. Otherwise typically about 50% of the mass would be lower than the mean, and 50% higher. (log(p) should be asymptotically normally distributed for large dimension d)
Here’s some Julia code to display what I mean. I’m pretty sure you’re aware of this math, but Anon and other readers may not understand it. Sometimes it’s just easier to explain with code and graphs:
Here’s the code in Julia:
https://gist.github.com/dlakelan/9949cdaf8c3b4bb67fe3d7e21e1cdb43
and here’s the resulting graph histogram/bar graph:
https://postimg.cc/cgz4sZqH
You can see that the “perfect” sample from a multivariate unit normal of dimension 20 has log(p) values that have a mean near about -28.0 and the lp of the mode is about -18.4 but the samples range from lp value about -40 to about -20
So if you’re at the mode, the point is “weird” because there’s very little volume that close to the mode and so the associated probability (density times volume) is small.
If you sampled uniformly from the distribution proportional to “indicator function on the set {x : lp(x) between -30 and -18.4} then you’d get a decent approximation to the posterior. If you then took those points and did a small number of diffusive metropolis steps to each one, you’d have essentially a perfect sample. If you did uniform on the region lp is -30 to -22 you’d to essentially just as well, basically there’s no real contribution from the region in the vicinity of lp = -18.4 because the region is small.
The “size of the region at a certain distance” is proportional to x^20, and I plot that at the end of the script, that’s like:
https://postimg.cc/9DvNGnGk
Volume grows extremely fast with distance. density declines rapdily with distance, the two effects combine to produce a region that has most of the mass that’s at lp values in a narrow range.
> Carlos, if you say “not lower than mean(log(p))-epsilon” for a sufficiently large epsilon then yes.
Right, “not lower than” is better rephrased as “not much lower than”. Just a little abuse of language – like in your example where -40 and -20 are “close” to -28 (but -18 is not).
> So if you’re at the mode, the point is “weird” because there’s very little volume that close to the mode and so the associated probability (density times volume) is small.
In that sense every other point is “weirder” because the associated probability is smaller.
> If you sampled uniformly from the distribution proportional to “indicator function on the set {x : lp(x) between -30 and -18.4} then you’d get a decent approximation to the posterior. […] If you did uniform on the region lp is -30 to -22 you’d to essentially just as well
Maybe I misunderstand what you wrote but I wouldn’t say that sampling “uniformly” gives a “decent approximation” when the probability density varies by several orders of magnitude within that region.
For {x: lp(x) between -40 and -22} the density is 65’659’969 times higher in some places than in others. (Looking at the histogram I suspect you wrote -30 when you meant -40, there is a relative factor in the thousands anyway.)
Carlos, no I did mean -30. I looked at the histogram of the log(p) and said “where is it high density?” and that was kind of -30 to -22 or something. If you add the region with even higher log(p) it doesn’t much matter because the region is very small, if you add the region where it’s well below -30 you’re adding way too much mass in the tail.
The thing is if you have a point around log(p) = -30 then a tiny euclidean perturbation in any direction is likely to drop log(p) by a bunch (there’s much more volume “farther away” from the mode than there is closer to the mode, so any perturbation at all is much more likely to take you “downhill” and the hill is likely to be steep). So in euclidean distance, only a small adjustment is needed to each sample point to equilibriate the sample.
So, sampling uniform on the region with log(p) greater than say -30, plus a few rounds of diffusive adjustments to each sample using metropolis hastings, and voila you have a very high quality sample.
The proof is in the pudding so to speak. Here’s the message on the stan forum from 2020 where I gave some output from my suggested sampler for a MvNormal with 100 dimensions.
https://discourse.mc-stan.org/t/the-typical-set-and-its-relevance-to-bayesian-computation/17174/22
The first set of graphs is for uniform on lp greater than quantile(lpvals, 0.99), which of course has light tails, and the second set is on lp greater than quantile(lpvals, 0.01) which of course has heavy tails.
When I think about what’s ideal, probably you want to sample light tails (not too light but perhaps a quantile of lp like .5 or .7 or .9), followed by diffusion, because diffusion will tend to spread things out and if you’re already too spread out it’ll require more luck to shrink things uphill, whereas if you’re too light tailed then each cycle of diffusion will spread you out in a way that likely improves the sample. It’s always better if your algorithm involving a random perturbation has a good probability of improving things.
On that basis perhaps you’d want to sample from uniform on lp greater than -22 or similar for my 20 dimensional example, and then diffusively equilibriate after that.
> So, sampling uniform on the region with log(p) greater than say -30, plus a few rounds of diffusive adjustments to each sample using metropolis hastings, and voila you have a very high quality sample.
That’s quite different from “a uniform distribution is a decent approximation of a normal distribution”.
(Also, I have the impression that the discussion in this post is precisely about how the Metropolis-Hastings algorithm has difficulties going “downhill” when the hill is steep. However I’ve not read all that has been written and I may miss something so I’ll leave it there.)
(It seems to me that in fact to get that uniform “decent approximation” to converge to the normal distribution you will need for those points to go mostly “uphill”. Such a movement is also problematic: even though it would be accepted it is unlikely to be proposed in the first place.)
Carlos, yes there are some subtleties. Hopefully I can code some stuff up in Julia and give some examples once I understand the API interface.
The proposal to use two types of samplers is motivated by the following thought:
1) It’s hard to move concertedly long distances using HMC when you are near the mode because you must take small steps due to the steep potential energy
2) If you are in a not terrible place which is fairly consistently too high log(p) then it’s easy to move small distances via diffusion to get into a more equilibriated location, since small changes to x will result in dropping the log(p) in most directions, and your limitation is only the size of the step to keep things transitioning. Yet, because of the steepness of the log(p) you probably don’t need large motions.
So the idea is start with 1) no gradient at all, just sample uniformly inside the appropriate boundary, moving in straight lines until you hit a boundary, potentially long distances, followed by 2) equilibriating this sample by diffusion with small step sizes to find your proper location on the fuzzy boundary.
When you try to do HMC you get potentially movement long distances, but to stay stable you need potentially small step sizes leading to a lot of function and gradient evaluations. However you only need the small step sizes in the presence of large gradients, which is not everywhere, just near the boundary… So you split the computation into two parts. One moving a long distance in straight lines with no computational difficulty and guaranteed acceptance, and one moving short distances in random directions. Neither stage has need for gradient calculations.
Metropolis-Hastings has a hard time going from high density to low density (or low energy to high energy, with potential energy being negative log likelihood in our fictitious physical systems). The reason is that the accept probability is min(1, p(proposal) / p(starting-point)). When the starting point is the mode, the proposals all have lower density.
I’ll try a different approach.
If you have the correct model, and happen to find the set of parameter values yielding the exact max density (I guess the “mode”), shouldn’t this be very informative about the “typical set”? Instead the algorithm is returning something that looks like it has zero information.
Say Omniscient Jones tells me in a dream the exactly correct model to predict where and when the next hurricane will form. Unfortunately it has 100k parameters that require fitting. If I run MCMC and it happens to come across the exactly correct set of parameter values (that allow us to know exactly where/when the next hurricane will form), it will either:
1) Blow up if we have chosen the optimal proposal distribution, which allows immediate convergence
2) Take forever due to very small step size
So essentially we cannot predict exactly where/when the next hurricane will form, even in principle, using this method.
Simulate some data points from a mixture of two gaussians with known parameters
mu_1=0.5, sigma_1=1
mu_2=-0.5, sigma_2=1
Let x_n be the value of exactly one simulated data point. Let’s try and fit this. Set our initial values
mu_1_inference = 0, sigma_1-inference = 5
mu_2_inference = x_n, sigma_2 inference = 1e-12
The density is extremely high. Indeed, it can be arbitrarily high as sigma_2_inference approaches zero. Do you think that makes it a perfect inference?
No, it wouldn’t. “The exactly correct set of parameters” would look approximately like all the other sets of parameters and almost certainly wouldn’t be the mode of your posterior. If you simulate 10 data points from normal(0, 1), then fit to it, the mode will not be (0, 1) because the data are insufficient to identify; you only know in this case that (0,1) is special based on information not in the data or model and which you wouldn’t have access to in a realistic scenario.
As the data gets larger, the posterior *mean* will converge in probability to identify the true value—not generally the mode. If you can prove for your posterior that mean = mode, which exists and is unique, and that your data are “large enough” that a single “best guess of the true parameters” suffices for your purposes, then yeah running MCMC will waste a bunch of your time. But the original post doesn’t say anything about those conditions.
Anon…
Remember that we have some finite dataset to work with, let’s say 1 hour of hurricane tracking data. Even if there is a perfect model from Omniscient Jones, and that perfect model has some perfect set of parameter values which would correctly predict everything… we have no way of knowing that based only on our 1 hour of hurricane tracking. Typically there will be a large set of parameters which will match the 1 hour of data all about equally well. If Omniscient Jones tells you to plug in a set of parameters but doesn’t tell you if they’re the perfect ones or just some random ones… can you tell which? No because lots of other parameter values will look about equally good based only on your 1 hour of data.
Now if OJ hands you a decade of minute-by-minute hurricane tracking data and an ultra-super-computer, you should ultimately wind up very very close to the perfect set of parameters that OJ tells you is optimal. But to do so will in fact require very very small step sizes because the posterior distribution will be extremely tightly constrained around the “perfect” values.
Thank you both. So then we are talking about overfitting. Ie, another condition is (roughly) that n_param/n_data ~ 0:
https://en.wikipedia.org/wiki/Double_descent
Possibly. I’m not sure what you mean. If OJ gives you the correct model and the correct parameter values, then by assumption, this is not over-fit but rather *correct*. But you are unlikely to learn those correct parameters from a realistic amount of data without OJ’s help.
On the other hand, if you gather some data and an incorrect but reasonable approximate model, and it just so happens that you have a small enough number of data points that your reasonable model can fit that data almost exactly… then yes you’re likely to be overfitting if you take just an optimal/modal set of parameters. However if you do an MCMC sample then you’re much more likely to get a range of values that shows you more realistic uncertainty.
That should have been “n_param/n_data ~ 1”.
But rather than assumption, we would know by correctly predicting data not used to fit the model.
Yes, this is very unlikely to occur. But if we do happen to find it (either during initialization/adaptation or fitting), it appears that MCMC will not tell you that you essentially won the lottery. Rather it returns something that looks like an error.
Seems to me instead you should be able to immediately jump to the “typical set” using this information.
Anon,
Note that if you have a lot of data, and you are concentrating about Omniscient Jones’ parameter values, that you’ll still have trouble initializing at the mode, but that the values you want to sample from are, in some sense, very close to the mode anyway.
One issue is how you measure distance in the parameter space. Let’s say that you have one parameter and its posterior is in essence normal(1.0,0.1) then you want to sample from values around 0.8 to 1.2
suppose instead you have normal(1.0,0.0001) then you want to sample from values around .9998 to 1.0002
If you measure in terms of “number of standard deviations” both are *the same* if you measure in terms of “fractions of the mean” (which is also just euclidean distance in parameter space in this case) then they’re very very different.
As you converge to a delta function you continue to have the issue of difficulty of sampling when you start at the mode, but how far you have to move drops dramatically the tighter the distribution.
If you try to sample normal(1.0,.0001) it requires very small steps, on the order of 0.0001, and most of the samples won’t be at 1.0 but they will be a very small distance away.
In some ways the same thing occurs in high dimensional space. You might just need to move not very far away (euclidean distance) from the mode to get into the “typical set” but it will result in a large change in probability density.
If you have already identified the correct model and “origin” (as shown by correctly predicting non-training data) why not move ~sqrt(d) away from it in some direction?
https://statmodeling.stat.columbia.edu/2020/08/02/the-typical-set-and-its-relevance-to-bayesian-computation/
And since you have already identified the target distribution, afaict now you also know the optimal sampling distribution.
The point is if a scientist does everything right (correct model, correct parameters, plenty of good data) the fitting method should not return what looks like an error.
Think of it like a unit test.
If it “looks like an error” it may because you still don’t understand what’s going on after this long exchange (that I’ve not read entirely). Or maybe because the distribution that you get represents indeed the “error” in taking a single value as answer when there is uncertainty.
> (correct model, correct parameters, plenty of good data)
I’m not sure what “correct parameters” means there, as I understand that the parameters are what you try to identify with the “fitting method” you mentions.
It’s not an error that the “fitting method” returns a distribution for the parameters of interest. It’s what you want when their value cannot be precisely determined. When problems are well defined and there is no uncertainty one usually solves equations instead of bothering with statistics.
If it is known that the distribution of parameters has the form of a d-dimensional unit normal one doesn’t need to bother with MCMC sampling either. There are more efficient methods of sampling.
Also, the d-dimensional unit normal example has nothing to do with the case where you have “plenty of good data”. If for example you want to determine position of something the posterior distribution may be a three-dimensional normal but the variance won’t be one: the distribution collapses to single point when you have “infinitely good data”.
> Most of our models don’t have modes
They do have modes – the problem comes from having too many of them.
No, most of our models do not have modes. For example, Andrew and Merlin’s model for the last election cycle for the Economist to predict Republican vote share. The model has no modes, not multiple modes. In MacKay’s book, he describes the problem as “EM goes boom.” The reason lme4 does max marginal likelihood rather than max likelihood is that hierarchical models don’t have modes.
Rather than modes, hierarchical models have poles as hierarchical variance goes to zero and the coefficients go toward the hierarchical mean (usually zero). Radford Neal’s funnel example is the purest—that’s a hierarchical model with no data.
Now you do run into problems with things like high-dimensional Gaussian mixtures (i.e., “soft” K-means) or neural network models. With simple priors/regularizers, those have multiple modes.
Bob, I always use zero avoiding priors for variances, usually gamma. I think of gamma as parameterized by a mode and a concentration parameter gamma(n,m/(n-1)). With m the mode and n the concentration. Typically concentration of 3 or 4 is my minimum, though I often use 10-20 depending on my knowledge.
This should prevent that funnel like behavior.
In a lot of my problems I’ve found the model fits really well and the problem is that the posterior for some parameters is really concentrated and therefore the step size must be very small.
> No, most of our models do not have modes.
I stand corrected, thanks for the clarification. Your models may not have modes. I was thinking about complicated likelihoods but didn’t consider divergences.
We can also have them in simple models – for example doing inference about a normal(mu, sigma) model from a single data point – and I see how things like adding variances can lead to unbounded likelihood functions in general.
I understand that when you write “hierarchical models don’t have modes” you are not categorically denying that they could. It seems that in principle one may use appropriate priors to tame the divergences – as Daniel suggested – or reparametrize the model to “fix” the likelihood (for example using precisions instead of variances).