Chain Ensembles

[Michael Betancourt]

Peter brought up the question of how to interface ensemble samplers with Stan, but if you step back the real question is how to interface ensembles of chains in Stan (in the simplest case an ensemble sampler is a single chain with all of other walkers just part of the proposal mechanism; it's easy to assign a single chain to each walker and then have them mix in the proposal stage, but this just generates more chains, not more samples).  So, just how do we want people to do inference with multiple chains?  

One option is to interleave the chains (carefully as to preserve detailed balance) and then do inference on the mega-chain.  My concern here is that the mega-chain might have a very weird autocorrelation structure (aliasing, single-frequency modes) that would invalidate most autocorrelation estimators (can't cut off the lags past the first zero anymore).

The other option is to calculate expectations from each chain, and then average them together using the ESS as weights.  This approach is more awkward for the user, as the averaging of the expectations need not commute with other manipulations (marginalization, etc).

Thoughts?

Once we've decided on this it won't be hard to modify command appropriately (for example, adding a chain index as a sampler parameter for ensemble methods or just having multiple output files).

[Andrew Gelman]

I'm not sure, but just in general I think it's important for our paradigm to be that we return to the user a set of simulations which he or she can use to represent the posterior distribution.  Calculating expectations is fine but that's just one thing that can be done with the simulations, I don't think it should be the paradigmatic activity.  And I really really don't want to be handing out weights--not unless there's a really really good reason to do so!

A

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Michael Betancourt]

Agreed -- giving samples to the user is the ideal.  But how can we aggregate --dependent-- samples?

[Andrew Gelman]

Regular MCMC samples are dependent and we aggregate them.

To put it another way, imagine the goal was to compute expectations.  In that case, any weighting of chains would be valid.  Equal weighting might not be the most efficient but it is still kosher.

One option is to do different rates of thinning of different chains.

A

[Michael Betancourt]

But is an aggregation of samples equivalent to a reversible Markov chain?  For example, if we

1) randomly sample a walker

2) evolve it 

3) take the evolved state as a sample

then we can consider the process as a mixture of Markov transitions.  Because each transition converges to the same stationary distribution, the mixture will, too.

But if we

1) evolve all walkers

2) append all evolved states as samples

then I don't think we have the same guarantee.  I'm not so much worried about convergence but the behavior of the autocorrelation estimators, in particular when/how we can terminate the lags in computing the ESS.

[Andrew Gelman]

Yes, I agree that for ESS we would like to keep the chains separate in the analyses.

Which reminds me that it would be cool if we could develop some sort of FFT version of ESS and split R-hat.

Which reminds me that we should write up our ESS and split R-hat stuff as a paper (currently it exists only in BDA3, the Stan manual, and in a 2-year-old draft paper by Ken Shirley).

A

[Bob Carpenter]

On 5/24/13 12:54 PM, Michael Betancourt wrote:
> But is an aggregation of samples equivalent to a reversible Markov chain?

I don't think so.

> For example, if we
>
> 1) randomly sample a walker
> 2) evolve it
> 3) take the evolved state as a sample
>
> then we can consider the process as a mixture of Markov transitions. Because each transition converges to the same
> stationary distribution, the mixture will, too.

Right.

> But if we
>
> 1) evolve all walkers
> 2) append all evolved states as samples
>
> then I don't think we have the same guarantee. I'm not so much worried about convergence but the behavior of the
> autocorrelation estimators, in particular when/how we can terminate the lags in computing the ESS.

I don't think we do, either.

For convergence and ESS, we want to keep the ensemble members separate.

For inference and posterior estimation, we want to throw them all together.
Andrew also wants to randomize the thrown-together samples to protect vs. client
code that might do something non-random, like take the first 200.

If the ensemble goes long enough, then the effect of just laying the
ensemble members end to end will diminish, so it should be the same in
the limit. But we hope to stop long before that happens! So what would
happen is that we'd underestimate autocorrelation by crossing ensemble
member boundaries.

- Bob

[Bob Carpenter]

On 5/24/13 12:57 PM, Andrew Gelman wrote:
> Yes, I agree that for ESS we would like to keep the chains separate in the analyses.
>
> Which reminds me that it would be cool if we could develop some sort of FFT version of ESS and split R-hat.

Something other than using an FFT in the implementation?
We're doing that already.

> Which reminds me that we should write up our ESS and split R-hat stuff as a paper (currently it exists only in BDA3, the
> Stan manual, and in a 2-year-old draft paper by Ken Shirley).

:-)

- Bob

[Andrew Gelman]

ESS is already OK that way, I guess.  But split R-hat is still a bit ugly, on a conceptual level.  I'd like something fully fractal if possible.

[Bob Carpenter]

> On May 24, 2013, at 1:10 PM, Bob Carpenter wrote:
>
>> On 5/24/13 12:57 PM, Andrew Gelman wrote:
>>> Yes, I agree that for ESS we would like to keep the chains separate in the analyses.
>>>
>>> Which reminds me that it would be cool if we could develop some sort of FFT version of ESS and split R-hat.
>>
>> Something other than using an FFT in the implementation?
>> We're doing that already.

On 5/24/13 1:11 PM, Andrew Gelman wrote:

> ESS is already OK that way, I guess.

Let us know if there's a better way to do this. Goodman has
an approach I haven't dove into in his ACOR code:

http://www.math.nyu.edu/faculty/goodman/software/acor/

which is also part of the Python package Dan Foreman-Mackey
used in the emcee sampler:

https://pypi.python.org/pypi/acor

> But split R-hat is still a bit ugly, on a conceptual level. I'd like something
> fully fractal if possible.

I see, so it's not how to compute split R-hat, but a new R-hat that
would involve splitting into 2, splitting into 3 (or 4), splitting
into 4 (or 8), and so on. No idea how to do it, but it does remind
me of how FFT works.

Too bad all our theory people are busy working on code. Your call is
important to us and we'll get to it as soon as our next agent is free :-)

(By the way, it's very odd getting e-mail from you before 4 PM --- I can
feel the time difference.)

- Bob

[Michael Betancourt]

I see, so it's not how to compute split R-hat, but a new R-hat that

would involve splitting into 2, splitting into 3 (or 4), splitting
into 4 (or 8), and so on.  No idea how to do it, but it does remind
me of how FFT works.

We can use DFT methods to compute relevant sums, but I don't think they're quite appropriate for Rhat.  The problem is that they give you they don't give you split partitions (first n/2 samples, last n/2 samples) but rather offset partitions (all even-indexed samples, all odd-indexed samples) so you don't get the same sensitivity to drift as you do in the split-Rhat case.  Unless I'm looking at it incorrectly?

[Bob Carpenter]

D'oh! You're right --- the splits in FFT work as you say, not
as in my weak analogy.

- Bob

[Michael Betancourt]

While waiting for acid_tests to run I spent some time looking through the ensemble literature, eventually going back to the original Goodman and Weare paper on which emcee is based.  Unlike the emcee paper (on the arXiv, glosses over waaaaaay too many details) Goodman and Weare treat the ensemble sampler as a Markov chain on the product space pi cross pi .... pi cross pi, and explicitly note that any the sequence of any particular walker is very much NOT a Markov chain.  In the end, though, they get tricky.  When comparing ESS between a single chain and the ensemble, they collapse the ensemble at each sample by averaging each walker to produce f_{i} = sum_{j} f (x_{ij}), where i indexes time and j the walkers.  They then compare ESS/accuracy using both 1D chains.

So one possibility for Stan is to output not output the full product space in samples.csv (or have multiple output files, one for each walker) and then have bin/print produce any summary stats using the above method.  Any nonlinear expectations needed by the user can be accommodated by adding transformed parameters in the original model.

Thoughts?

- Bob

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Bob Carpenter]

On 6/13/13 12:54 PM, Michael Betancourt wrote:
> While waiting for acid_tests to run I spent some time looking through the ensemble literature, eventually going back to
> the original Goodman and Weare paper on which emcee is based. Unlike the emcee paper (on the arXiv, glosses over
> waaaaaay too many details) Goodman and Weare treat the ensemble sampler as a Markov chain on the product space pi cross
> pi .... pi cross pi, and explicitly note that any the sequence of any particular walker is very much NOT a Markov chain.

Right. Their equation (5) states (but does not explain why)
that for J walkers at iteration i, you have:

p(x_{i,1},...,x{i,J}) = p(x_{i,1}) ... p(x_{i,J})

I'd like to understand why. It seems wrong to me given
that x_{i,j} depends on x_{i-1,j} and also at least two
other points, other x_{i-1,j'} and x_{i-1,j''}.

> In the end, though, they get tricky. When comparing ESS between a single chain and the ensemble, they collapse the
> ensemble at each sample by averaging each walker to produce f_{i} = sum_{j} f (x_{ij}), where i indexes time and j the
> walkers. They then compare ESS/accuracy using both 1D chains.

In section 3, they seem to be restricting attention to estimating
posterior expectations means.

> So one possibility for Stan is to output not output the full product space in samples.csv (or have multiple output
> files, one for each walker) and then have bin/print produce any summary stats using the above method. Any nonlinear
> expectations needed by the user can be accommodated by adding transformed parameters in the original model.
> Thoughts?

We want to do full Bayesian inference, so just averaging over
walkers j isn't going to help us because the posterior
variance will be too low. But, if that can compute a reasonable
effective sample size that we could use with all the samples,
that could work.

But if each walker itself isn't a Markov chain, then we can't
even just use one chain. So I'm not sure what we can do with
the "samples."

- Bob

[Bob Carpenter]

OK -- caught that one. It's CTRL-n that sends the
friggin' message. I really wish people used the built-in
OS-specific functionality!



- Bob

> --
> You received this message because you are subscribed to the Google
> Groups "stan development mailing list" group.

> To unsubscribe from this group and stop receiving emails from it,
> send an email to stan-dev+u...@googlegroups.com.

[Michael Betancourt]

> Right. Their equation (5) states (but does not explain why)
> that for J walkers at iteration i, you have:
>
> p(x_{i,1},...,x{i,J}) = p(x_{i,1}) ... p(x_{i,J})
>
> I'd like to understand why. It seems wrong to me given
> that x_{i,j} depends on x_{i-1,j} and also at least two
> other points, other x_{i-1,j'} and x_{i-1,j''}.

Yup, each point in the ensemble depends on each other -- that's why the individual sequences are not Markov chains. But if you treat the ensemble as a state, {x_{0, j}, …, x_{i, j}, … x_{I, j}} then the distribution of the ensemble is preserved by each move. Basically you're parameters are the parameters of the full state and you're just Gibbs sampling when you update one of the walkers.

> In section 3, they seem to be restricting attention to estimating
> posterior expectations means.

Which, to be fair, is about all you can do with MCMC samples.

> We want to do full Bayesian inference, so just averaging over
> walkers j isn't going to help us because the posterior
> variance will be too low. But, if that can compute a reasonable
> effective sample size that we could use with all the samples,
> that could work.

But you don't use the posterior variance of that average. You average x_{j} = 1 / I \sum_{i} x_{ij}, and the sequence { x_{j} } becomes a valid Markov chain. You compute your expectations, variances, ESS, and MCSE from those averages. In fact it's actually more conservative than treating each chain as independent and weighting the samples by the individual chain ESSs.

> But if each walker itself isn't a Markov chain, then we can't
> even just use one chain. So I'm not sure what we can do with
> the "samples."

If you just took one walker and threw everything else out, then you'd have a Markov chain. If you repeated the same process with each walker you'd have something that looked like a Markov chain but together they'd be correlated and so you couldn't treat them independently. By averaging you get around the correlations, but become sensitive to chains that haven't converged (I think you can get around this by using the median).

Basically, you can't use each individual walker as a separate draw. But if you map ( x_{0, j}, …, x_{i, j}, … x_{I, j} ) -> x_{j}, say via an average or a median, for a given time index j then you get a sequence of draws that is a proper Markov chain.

[Bob Carpenter]

Why is the average of the J walkers a Markov chain? If
mu[i] is the average at i, I can't see how to formulate

p(mu[i] | mu[i-1])

Is there something general that says a transform of a Markov
chain is a Markov chain? If so, that'd explain why both the
average and the projection to one walker form Markov chains.

I don't see how you can estimate posterior variance from
the mu[i], though. It seems to me that the variance of
the chain of mu[i] will be lower than that of the component
walkers.

- Bob


P.S. Thanks for trying to explain this all to me.

[Michael Betancourt]

> Why is the average of the J walkers a Markov chain? If
> mu[i] is the average at i, I can't see how to formulate
>
> p(mu[i] | mu[i-1])
>
> Is there something general that says a transform of a Markov
> chain is a Markov chain? If so, that'd explain why both the
> average and the projection to one walker form Markov chains.

Just consider the multivariate case: you have a Markov chain sampling from some p(x1, x2). You can transform to p( x2 ) = int dx1 p(x1, x2) at each state, so that { x2 } are valid draws from a Markov chain for p( x2 ).

The proof is straightforward. Let T be a Markov transition with invariant density pi(x1, x2). Then by definition

pi(x1, x2) = int dx1' dx2' T (x1, x2 | x1', x2' ) pi(x1', x2').

Integrating both sides with respect to x2 gives

int dx2 pi(x1, x2) = int dx2 dx1' dx2' T (x1, x2 | x1', x2' ) pi(x1', x2')

or

pi(x1) = int dx1' [ int dx2 dx2' T(x1, x2 | x1', x2') pi(x1', x2') / pi(x1') ] pi(x1')

pi(x1) = int dx1' [ int dx2 dx2' T(x1, x2 | x1', x2') pi(x2' | x1') ] pi(x1')

Where

T(x1 | x1') \equiv int dx2 dx2' T(x1, x2 | x1', x2') pi(x2' | x1')

is a valid Markov transition for pi(x1).

Neal talks about the general case in his ensemble paper, and the intuition is really nice. The marginalized transition first samples x2 from the condition distribution p(x2 | x1) and then performs a transition on the space (x1, x2) before marginalizing back to just x1. It's also how HMC works, where x2 is the momenta, x1 the position, and T(x1, x2 | x1', x2') the Hamiltonian flow + Metropolis correction.

You can readily generalize the proof to include changes of variables x1, x2 -> f(x1, x2) (where f(x1, x2) might be the average or the median), and the same reasoning holds for the ensemble where you consider the state space x = {x1, …, xi, … xI}.

> I don't see how you can estimate posterior variance from
> the mu[i], though. It seems to me that the variance of
> the chain of mu[i] will be lower than that of the component
> walkers.

Right, but you only use the mu[i] to compute the expectation of mu. If you want the variance then you have to compute var[i, j] in each walker and average that. While this method does account for the correlations in a systematic way it is terribly inflexible.

> P.S. Thanks for trying to explain this all to me.

No worries at all. If people didn't push me to explain I'd get in the habit of handwaving instead of whipping up little proofs like that and understanding it better myself. Keep me on my toes!

[Michael Betancourt]

>> I don't see how you can estimate posterior variance from
>> the mu[i], though. It seems to me that the variance of
>> the chain of mu[i] will be lower than that of the component
>> walkers.
>
> Right, but you only use the mu[i] to compute the expectation of mu. If you want the variance then you have to compute var[i, j] in each walker and average that. While this method does account for the correlations in a systematic way it is terribly inflexible.

Actually, this is exactly the behavior we want! The variance of the means will indeed shrink, which means the MCSE will also shrink. But that is the desired behavior -- if the walkers are independent then the MCSE will decrease by 1 / n_walkers, giving MCSE = var / (ESS * n_walkers). The full power of the samples is preserved.

[Michael Betancourt]

I TeX'd the derivation up as it was a big ungainly in plain text.

[Bob Carpenter]

On June 14, 2013 06:48:37 AM EDT, Michael Betancourt <betan...@gmail.com> wrote:

> I TeX'd the derivation up as it was a big ungainly in plain text. Thanks!

Thanks. That makes things much clearer. As to the note,
I think it might help to discuss the estimation of variance
as an example in the last paragraph. But don't bother on my
behalf --- I see what's going on now.

It also made me realize that the n_eff we report
is the n_eff for the mean, so that won't carry over to estimating
variance, either.

Could we just use the ensemble averages to compute n_eff, then

just go ahead and use all the samples for inference? Each of the samples
has the right marginal distribution. We won't necessarily
know the n_eff, but then we don't know that for what we have anyway
unless we compute it on the outside.

- Bob

[Michael Betancourt]

> Thanks. That makes things much clearer. As to the note,
> I think it might help to discuss the estimation of variance
> as an example in the last paragraph. But don't bother on my
> behalf --- I see what's going on now.

Yeah, it would be an important for any implementation.

> It also made me realize that the n_eff we report
> is the n_eff for the mean, so that won't carry over to estimating
> variance, either.

Nope, for that you'd need fourth order statistics which no one ever quotes. I've never seen more than MCSE quoted or used.

> Could we just use the ensemble averages to compute n_eff, then
>
> just go ahead and use all the samples for inference? Each of the samples
> has the right marginal distribution. We won't necessarily
> know the n_eff, but then we don't know that for what we have anyway
> unless we compute it on the outside.

You could, but then why did we calculate n_eff? Correlations don't effect the mean, just the variance/MCSE so users could use ensemble samples or ensemble averages to compute whatever expectation they desire, they'd just have no good idea of its precision. If we're trying to push users into more "responsible" use of MCMC then we'd have to force them into using that they specified and was averaged, if it's not defined in your model or in the transformed parameters block then we can't make any claim of confidence and so won't let you use it. At least not use it easily.

[Andrew Gelman]

I see n_eff as having three roles:
1. We use n_eff in our (informal) stopping rules, for example we require n_eff > 10 for all quantities of interest.
2. n_eff is used in computing the Monte Carlo standard errors
3. We use n_eff in our efficiency comparisons (seconds of wall time per n_eff).

[Bob Carpenter]

On 6/14/13 10:21 PM, Andrew Gelman wrote:> I see n_eff as having three roles:

This is exactly what I'm worried about. And not
just for ensemble methods.

We report posterior sd and posterior quantiles, but
we only report n_eff for the mean. What's our
certainty on the posterior quantile and sd estimates?

We can calculate theta^2 in a single sample, but not
variance or sd or quantiles, because they depend on
multiple samples. In comparing systems, Matt had
been using theta^2 for all variables theta of interest
to get a handle on how well posterior variance was being
estimated.

Other expectations, like event probabilities, can be pretty
easily calculated by defining variables using indicator
functions. These should have their uncertainties calculated
properly.

- Bob

[Michael Betancourt]

> > 1. We use n_eff in our (informal) stopping rules, for example we require n_eff > 10 for all quantities of interest.
> > 2. n_eff is used in computing the Monte Carlo standard errors
> > 3. We use n_eff in our efficiency comparisons (seconds of wall time per n_eff).

Agreed. The only issue is that all three require knowing the expectation of interest a priori.

> This is exactly what I'm worried about. And not
> just for ensemble methods.
>
> We report posterior sd and posterior quantiles, but
> we only report n_eff for the mean. What's our
> certainty on the posterior quantile and sd estimates?

I'm not sure about the quantiles, but you can compute the variance of the variance and define a MCSE the same as for the mean. In practice, however, the uncertainty o the variance isn't as critical -- it's the uncertainty of the posterior uncertainty and second-order in some sense.

> We can calculate theta^2 in a single sample, but not
> variance or sd or quantiles, because they depend on
> multiple samples. In comparing systems, Matt had
> been using theta^2 for all variables theta of interest
> to get a handle on how well posterior variance was being
> estimated.

Right -- technically the variance is a function of the data and an unknown constant (the mean). It can be calculated via the expectations of theta and theta^{2}, however, which would be straightforward to implement in both the single state and ensemble schemes.

> Other expectations, like event probabilities, can be pretty
> easily calculated by defining variables using indicator
> functions. These should have their uncertainties calculated
> properly.

Here's the way I see the grand scheme of things. One of the great features of MCMC is that the samples can be transformed or marginalized trivially, without changing their validity. That means you can run your sampler and still have valid inference even if you decide that you're interested in different expectations. On the other hand, this flexibility requires the user to compute ESS themselves. We can handle the calculations automatically, provided the user define the desired expectations a priori. This leave three approaches

One Stage Stan:

All expectations have to be defined in the model or transformed parameters. We run the model, compute mean, variance, ESS, and MCSE automatically and report those summary statistics.

One Point Five Stage Stan:

We compute the summary statistics automatically as above, but store the intermediate samples in case an expert user wants to change the expectations.

Two Stage Stan:

Sampling and post-processing are two completely different stages, and we rewrite bin/print to be less reliant on the Stan format and assume little more than a csv.

Currently, Stan uses the middle strategy. One Stage Stan would be more user friendly but slightly more annoying for expert users. Of course, it wouldn't be hard to make the approach customizable.

[Ben Goodrich]

On Saturday, June 15, 2013 11:00:38 AM UTC-4, Bob Carpenter wrote:

This is exactly what I'm worried about.  And not
just for ensemble methods.

We report posterior sd and posterior quantiles, but
we only report n_eff for the mean.  What's our
certainty on the posterior quantile and sd estimates?

And the n_eff is only applicable to the marginal mean. Perhaps that is okay if the user were only interested in some treatment effect. But it really doesn't make sense for the multivariate parameters such as the simplex, correlation matrix, etc. where a user pretty much can't be interested in the marginal means in isolation. It would be nice if we had some measure of effective sample size for a multivariate parameter as a whole, but in the absence of that, I guess we should be telling users to calculate all quantities of interest in the generated quantities block.

Ben

[Andrew Gelman]

On Jun 15, 2013, at 11:00 AM, Bob Carpenter wrote:

We report posterior sd and posterior quantiles, but
we only report n_eff for the mean.  What's our
certainty on the posterior quantile and sd estimates?

I don't see the need, in most cases, to give uncertainties on posterior quantiles and sd's.  If these are desired, we can do so (most simply by computing them separately for each of the split chains and then getting uncertainties that way, more generally by some sort of fft-like approach).  But most of the time I think the important things are the parameters themselves, not the quantiles etc.