parameter chains mixing well in some regions, stuck in others

[Jeffrey Arnold]

This is probably due to a bad model, but I am getting some behavior in my chains that is sufficiently unusual that I am stumped. The overview is that the chain will switch between regions with good mixing behavior and regions in which the parameters do not move. The regions with good mixing have a high treedepth, while the regions in which the parameters are approximately constant have a low treedepth (2,3). 

I attached plots of the chain of the lp__ (all other parameters have qualitatively the same pattern) and the treedepth. These 1000 iterations were from a run of 40,000 with a warmup of 20,000.   I've been getting this behavior in various dynamic linear models that I've been running lately.   Since it is of non-trivial size, I put the the model, output, and input of a model that produced this behavior in a zip file here: https://dl.dropbox.com/u/482019/m_prices_1.zip.  

Sorry that I couldn't make it simpler, I don't expect you to work through the model; I was just wondering if you have any ideas of what sort of problems this chain behavior might indicate. 

-Jeff 

[Bob Carpenter]

Low tree depth and bad mixing probably
means the step size is too large in that region
and there's too much rejection going on.

Is this why you were asking about input of varying
step sizes? What I'd suggest doing is using smaller
step size(s), but as you pointed out in the last message,
we don't have a way to let you do that!

One thing you might try to get out of it is turn
on the epsilon_pm parameter and let the step size
get significantly smaller. It'll slow everything
else down, but it may allow you to mix in bad regions.

We've also been considering automatically lowering
the step size every 10th iteration or so.

Another thing that can help is reparameterizing so
the posterior space is better behaved. What we've
been calling the "Matt trick", because he introduced
it. There's an example in the manual for the new release.
But the basic idea in the simplest case
is that rather than doing this:

for (n in 1:N)
alpha[n] ~ normal(mu,sigma);

... use alpha[n] ...

do this:

for (n in 1:N)
alpha_raw[n] ~ normal(0,1);

... use (mu + sigma * alpha_raw[n]) ...

I really wish we could do this kind of transform automatically!

- Bob

> --
>
>

[Ben Goodrich]

On Tuesday, November 27, 2012 11:00:52 PM UTC-5, Jeffrey Arnold wrote:

Sorry that I couldn't make it simpler, I don't expect you to work through the model; I was just wondering if you have any ideas of what sort of problems this chain behavior might indicate. 

Bob stated the Matt trick in general, and there are a few examples of it on the list or in the manual. If you search for "random walk" on stan-users, there are a couple of things that are more pertinent to your situation. In particular, I think you would want to make a theta_innovation parameter with a standard normal prior and then make theta a transformed parameter that is theta[i] <- theta[i-1] + sigma * theta_innovation[i].

Also, the constraint that the loadings matrix is unit lower-triangular is convenient, but it does make all the free parameters sensitive to the first n_factors outcome variables. If those outcomes are not that germane, or have a small amount of common variance, or have very different variances, then you could get computational problems. So, you might try reordering the outcomes to put the theoretically best ones first.

Ben

[Jeffrey Arnold]

On Tuesday, November 27, 2012 11:21:20 PM UTC-5, Bob Carpenter wrote:

Low tree depth and bad mixing probably
means the step size is too large in that region
and there's too much rejection going on.

Is this why you were asking about input of varying
step sizes?  What I'd suggest doing is using smaller
step size(s), but as you pointed out in the last message,
we don't have a way to let you do that!

It was partially that :)

 

One thing you might try to get out of it is turn
on the epsilon_pm parameter and let the step size
get significantly smaller.  It'll slow everything
else down, but it may allow you to mix in bad regions.

When I set --epsilon=whatever,  is the result is that there is no adaption of step size or tuning of parameter specific step sizes, but the NUTS algorithm still works? And if there is no adaption, does anything happen in warmup other than simply discarding the iterations? 

 

Another thing that can help is reparameterizing so
the posterior space is better behaved.  What we've
been calling the "Matt trick", because he introduced
it.  There's an example in the manual for the new release.
But the basic idea in the simplest case
is that rather than doing this:

I went and read the things on the "Matt trick" and I'm definitely going to give it a try. 

If you have time, I was hoping for some more intuition why the Matt trick works? 

I thought that Stan was good with high correlation but bad with different scales of variables, which is why 

you added the parameter specific step sizes in NUTS II? Is it the combination of differing scales + high correlation  in parameters

that still is difficult for STAN? 

[Ben Goodrich]

On Wednesday, November 28, 2012 3:53:17 PM UTC-5, Jeffrey Arnold wrote:

Another thing that can help is reparameterizing so
the posterior space is better behaved.  What we've
been calling the "Matt trick", because he introduced
it.  There's an example in the manual for the new release.
But the basic idea in the simplest case
is that rather than doing this:

I went and read the things on the "Matt trick" and I'm definitely going to give it a try. 

If you have time, I was hoping for some more intuition why the Matt trick works? 

I thought that Stan was good with high correlation but bad with different scales of variables, which is why 

you added the parameter specific step sizes in NUTS II? Is it the combination of differing scales + high correlation  in parameters

that still is difficult for STAN? 

 In short yes, high correlations can exacerbate the problem of differing scales. I think Marcus gave the best explanation at

https://groups.google.com/d/msg/stan-users/2YMo3SUzHTg/4xfu1IKVT78J

but be sure to unhide the quoted text to see the rest of it.

If theta were a parameter and

theta[i] ~ normal(theta[i-1], sigma)

then the joint distribution of theta[i], theta[i-1], and sigma is very correlated. Conversely, if theta_innovation is parameter with iid standard normal priors, then they've got the same scale and no dependence, at least a priori. Then, theta is a transformed parameter such that

theta[i] <- theta[i-1] * sigma * theta_innovation[i]

which may affect the likelihood of the data but doesn't get a prior.

I think there is also a sense in which Stan is optimized for sampling from independent standard normals with the way the momentum and leapfrogging works, so whatever you can do to make the posterior distribution be closer to N(0,I) is good.

Ben

[Bob Carpenter]

It's not tied to normality or correlation per se.

If the eigenvalues are all the same, Stan will work equally
well no matter what the rotation is.

Basic HMC is rotation invariant, but it's not scale invariant.

Stan's varying step sizes rescale each dimension
which destroys the rotation invariance.

But the really big issue is the position-dependent geometry of
the posterior. Specifically, in a hierarchical model,
when the hierarchical variance is low, the scale of
the lower level parameters will be low, and when it's
high, the scale will be high. This is clearly illustrated
with Radford Neal's funnel density. I did a plot and provide
the reparameterization. I put the draft of the relevant section
up on out code page:

http://code.google.com/p/stan/downloads/detail?name=stan-reference-1.0.4-draft-reparam.pdf

As an added bonus, you get a section from Ben on multivariate
reparameterizations.

There's more to come in this section. The normal isn't the
only thing that can get reparameterized this way.

And, as Ben points out, we'd ideally like to standardize
the posterior, but we don't know what that is, or we
wouldn't be doing MCMC in the first place!

- Bob

[Bob Carpenter]

On 11/28/12 3:53 PM, Jeffrey Arnold wrote:
> On Tuesday, November 27, 2012 11:21:20 PM UTC-5, Bob Carpenter wrote:

...

> When I set --epsilon=whatever, is the result is that there is no adaption of step size or tuning of parameter
> specific step sizes, but the NUTS algorithm still works? And if there is no adaption, does anything happen in warmup
> other than simply discarding the iterations?

...

--epsilon=<double> pins the basic step size to that value.

--equal_step_sizes turns off varying step size adaptation.

Setting both of these will use --epsilon for the step size
in each dimension with no adaptation whatsoever.

--L=<int> sets the number of leapfrog steps (and turns off NUTS)

If you don't use --L, then NUTS determines the number of steps.

- Bob

[Jeffrey Arnold]

Thanks Ben and Bob. I have a somewhat better intuition now.  That part of the reference manual looks really good. 

This is somewhat related, but does parameter expansion have any effect in Stan? More generally, a list of which heuristics / tricks from gibbs sampling carry over to HMC would be very useful. 

[Bob Carpenter]

On 11/28/12 11:44 PM, Jeffrey Arnold wrote:
> Thanks Ben and Bob. I have a somewhat better intuition now. That part of the reference manual looks really good.
>
> This is somewhat related, but does parameter expansion have any effect in Stan?

I'll quote from Matt, who's much more articulate on
these matters:

From Matt Hoffman, 14 Aug 2012
=========================================
All I did was replace the hierarchical prior

σ ~ U(0, 1000)
μ ~ N(0, 1000)
θ_j ~ N(μ, σ)

with the equivalent non-hierarchical prior

σ ~ U(0, 1000)
μ ~ N(0, 1000)
\tilde θ_j ~ N(0, 1)
θ_j = \tilde θ_j * σ + μ.

Basically, it's just a change of variables that makes the nasty
heavy-tailed hierarchical prior into a simple product prior on a
different measure space. In graphical model terms, this reads as
taking apart the hierarchy and making all of the nodes point directly
at the data.

Empirically, the model as originally written fails to get into the
left tail of σ (and rejects a lot when it gets close), whereas the
reparameterized model has no such problems. Mixing also seems faster,
which isn't too surprising, given that the reparameterized model has a
much smaller range of log-probabilities (std. dev. of about 2.6)
compared with the original model (std. dev. of about 6.5).
=========================================

It's like the funnel only a more realistic example and Matt
throws in an explanation of why it works.

> More generally, a list of which
> heuristics / tricks from gibbs sampling carry over to HMC would be very useful.

We're still working a lot of that out ourselves!

What will help most for Stan is any transformation that
makes the posterior more like N(0,I).

Either scaling or rotation on their own should be OK,
but together they can be problematic.

In general, it helps to apply transforms if they

-- eliminate varying scales in different parts of the posterior
that would ideally be suited to different step sizes
-- anything that gets rid of correlation related to the above

The funnel transform is an example of both. Also,

-- anything that makes the scale more uniform
at varying points in the space

For the artificial funnel example, there's strong correlation
driving position-dependent scaling, which is easy to convert
to unit normal.

-- anything that makes thin tails or thick tails more normal

The above example from Matt is an example of thinning the tails,
and the triangle distribution (in the current v.1.0.3 manual)
by transforming through constraints -- that's an implicit logit
transform that changes the tail behavior.

Anyone else have anything else for the list? I can include
it in the reparameterization section of the manual.

- Bob

[Jeffrey Arnold]

The reparameterization has much better mixing. Thanks!

I've been doing a lot of dynamic linear model-type models in Stan. I was wondering if for DLMs, would it (a) make sense and (b) be feasible to have a Kalman filter function in Stan? 

For a set of latent state parameters, would it be more efficient to calculate the log posterior contribution using Kalman filter "distribution"? 

[Ben Goodrich]

On Friday, November 30, 2012 5:02:00 PM UTC-5, Jeffrey Arnold wrote:

The reparameterization has much better mixing. Thanks!

I've been doing a lot of dynamic linear model-type models in Stan. I was wondering if for DLMs, would it (a) make sense and (b) be feasible to have a Kalman filter function in Stan? 

For a set of latent state parameters, would it be more efficient to calculate the log posterior contribution using Kalman filter "distribution"? 

Do you have a specific link for your proposal? The first thing that comes to my mind are that the multivariate normal distribution is very slow in Stan, so it might be faster in theory but slower in practice (due to Stan's need to auto-diff it). Also, Stan would need something that pertains to a joint distribution rather than a conditional distribution of the latents given the observed variables.

Ben

[Jeffrey Arnold]

Alas, I don't have any answers or ideas at the moment. That question was more fishing for answers than being able to provide one. I suppose that something like a Kalman filter could be used to update the log posterior by marginalizing over the latents rather than explicitly calculating them. But I'd need to think about it more. 

--
 
 

[Ben Goodrich]

On Saturday, December 1, 2012 1:34:48 PM UTC-5, Jeffrey Arnold wrote:

Alas, I don't have any answers or ideas at the moment. That question was more fishing for answers than being able to provide one. I suppose that something like a Kalman filter could be used to update the log posterior by marginalizing over the latents rather than explicitly calculating them. But I'd need to think about it more. 

Maybe something like this?

http://scholar.google.com/scholar?cluster=17139112289009455747&hl=en&as_sdt=0,33

Ben

[Matt Hoffman]

Another option to lower step sizes is to use a larger value of --delta
(i.e., closer to 1). This will cause the step size adaptation to shoot
for more accurate simulations and therefore smaller step sizes. Should
be easier (and better calibrated) than trying to set epsilon by hand.

I changed the command-line doc to be a little clearer about what --delta does.

> --
>
>