A question on convergence

[Mike Boedigheimer]

Hi

I looked and I haven't seen the following result mentioned anywhere> I know you put up with a lot and I apologize if it is trodden ground. Two questions below will be is the result known? Is it something to worry about? 

I have fit a model with around 100 parameters on 8 threads with 1000 simulations per thread. Other than threads, I'm using all default parameters. R hat is 1 for all of the parameters and mixing looks good. No problems there.

The interesting part is if someone were to show me a "new" posterior draw from a one of these threads, I would be to predict which thread it came from with 60% accuracy (uniform positive predictive value across thread). That is pretty good compared to a random guess of 12.5%. I base this claim on 10-fold cross-validation using a quadratic discriminant classifier using the 100 parameters as predictors for the threads.  KNN also performs well.

The result persists using 1000 to 2000 or 4000 warmups. 

The result does not seem to be specific to any one dataset. But, it could be model specific. I've only tested one model which is an hierarchical binomial with random beta variables for within and between subject variability. I posted that model on an earlier post.

This result seems to imply that the samplers in each thread are exploring different parts of the posterior or are dwelling too long on particular spots. 

Honestly, the posterior means don't seem to change too much due to this, but they seem to change a little. If I compare the posterior means from one of these threads to  a "true" mean estimated from very sparse sampling.  

Thanks for all your help. 

[Ben Goodrich]

On Tuesday, October 27, 2015 at 8:26:58 PM UTC-4, Mike Boedigheimer wrote:

The result does not seem to be specific to any one dataset. But, it could be model specific. I've only tested one model which is an hierarchical binomial with random beta variables for within and between subject variability. I posted that model on an earlier post.

This result seems to imply that the samplers in each thread are exploring different parts of the posterior or are dwelling too long on particular spots. 

Honestly, the posterior means don't seem to change too much due to this, but they seem to change a little. If I compare the posterior means from one of these threads to  a "true" mean estimated from very sparse sampling.  

I'm not too surprised by this because the Rhat is really only telling you about the posterior means. Can you attach the whole .stan program and some real or simulated data?

Ben

[Mike Boedigheimer]

Sure. Attached are is the stan program and some data. Let me know what you find. Thanks! 

[Ben Goodrich]

On Wednesday, October 28, 2015 at 11:11:25 AM UTC-4, Mike Boedigheimer wrote:

Sure. Attached are is the stan program and some data. Let me know what you find. Thanks! 

I did 16 chains and they look fine collectively. What did you do to predict which chain a draw came from?

Ben

[Mike Boedigheimer]

That is what is so surprising. To make the predictions I used 10-fold crossvalidation. One each training set I fit a quadratic discriminant (k-nearest neighbors also works). With the full specified model I predicted which thread each of the observations in the out-sample came from. I repeated this for all 10 folds. If you have matlab I can share all that code.

I also looked at PCA of the (centered and unit-variance normalized) parameters and used permutations to assess the number of non-random components ( uses the method of buja and eyuboglu described in Leek and Story (PNAS, 2008 v105 p18718). There were more than 8. I think there should be 0 or 1 ( I need to think about that some more). 

I believe HMC is behaving like an old-school pseudo-random generator, 1st order is right on the mark, but there is some pattern in there. 

[Bob Carpenter]

That's expected for MCMC --- the draws in each Markov chain
are correlated. That's why the effective sample size
is less than the sample size.

If the 10-fold x-validation used a random subset of
points in each chain rather than removing a contiguous
10% subset, I wouldn't be surprised if the accuracy was
even higher from classification. I'd also guess that accuracy
goes down as the chains get longer. Imagine Metropolis
here --- you get rejections and then have copies of
exactly the same point in a chain. The probability of
that same duplication in another chain is measure 0.

- Bob

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

[Mike Boedigheimer]

Bob, I believe your predictions are true. I.e. if I train on the 1st half of each chain and predict the second half, I don' t get strong predictions (sometimes anti-predictive). It is also true that longer chains are slightly weaker. Although it is still easily detectable up to 40,000 total simulations. Also, if I thin the data by 1/8 the predictive value is not detectable. So I have to agree that correlation is a culprit. 

I will add one more thought by suggesting PCA might be a more sensitive way to detect this. I still find non-random latent variables in the PCA analysis after thinning the data by 1/40 or more. 

I am concluding that thinning is still best practice for HMC (even though I read many times that thinning should not be necessary with HMC). It may not matter much, or at all, for inference on a single parameter, but it looks like it might when making inference on linear combinations of the parameters (as I am in this problem). If you disagree please let me know. 

As always I want to thank you for your and the teams dedicated work on STAN and this forum.

[Ben Goodrich]

On Wednesday, October 28, 2015 at 7:44:03 PM UTC-4, Mike Boedigheimer wrote:

I am concluding that thinning is still best practice for HMC (even though I read many times that thinning should not be necessary with HMC). It may not matter much, or at all, for inference on a single parameter, but it looks like it might when making inference on linear combinations of the parameters (as I am in this problem). If you disagree please let me know. 

Thinning is not good for substantive inference because the draws contain some information even if they are not quite independent within chain. Your example seemed to mix really well though, so there was not very much within-chain dependence at all. Once in a great while (with NUTS) thinning is a necessary evil to work around hardware limitations. But the best thing to do is make the linear combination of stuff in generated quantities so that you can verify its MCMC standard error is acceptably small to make inferences with.

Ben

[Mike Boedigheimer]

Ok. I thought I was good, and now I'm confused. Are you saying that it is worse to thin 1e6 draws by drawing every 1/1000 samples than it is to accept chains with known dependence from 1000 draws?  how about 1 sample each from 1e3 chains? My intuition is way off here, so if you can point me to some reference materials that speak to this I'd really appreciate it. 

Even if I generate the linear combinations as quantities within STAN, I'm not sure I'm safe to trust MC standard error. If the draws are correlated it cannot be safe unless the chains are long enough that that the tails of all the interesting quantities are evenly sampled!

So, for now, I've updated my my position to think I may not have sampled enough so that the linear combinations of interest are evenly sampled. I will base my assessment on whether the PCA analysis reveals NO latent variables among the normalized parameters. I will see whether I can investigate what happens if I run 1000 threads with 1e9 samples each. Let me know if you think this is interesting or a (predictable) waste of time (and what the outcome will be;). Will the PCA analysis converge to 0 latent variables?

Best regards

MIke

[Andrew Gelman]

1.  If your sampler is moving so slowly that you would be throwing away 99.9% of your draws, you have problems.  It’s possible for this to be ok, but usually when we do HMC and it’s working, the efficiency is much higher than 1/1000.

2.  The effective sample size takes into account correlation of the draws.  You can use the effective sample size from Stan, it is not necessary for draws to be independent.

[Mike Boedigheimer]

Actually, in this case already I already have 100% efficiency (8000 Reff/8000 samples) for most parameters, Rhat is 1 and there is excellent mixing on a univariate level. I'm talking about eliminating any correlations by thinning after the simulation is done to remove thread-specific multi-parameter co-variance or means

[Ben Goodrich]

On Thursday, October 29, 2015 at 9:00:19 AM UTC-4, Mike Boedigheimer wrote:

Actually, in this case already I already have 100% efficiency (8000 Reff/8000 samples) for most parameters, Rhat is 1 and there is excellent mixing on a univariate level. I'm talking about eliminating any correlations by thinning after the simulation is done to remove thread-specific multi-parameter co-variance or means

Thinning is not going to help. When the output comes up and says that you have this or that ESS for a parameter, it indicates that the average of all the somewhat dependent draws has the same precision as you would expect from a simple random sample from the posterior distribution of size ESS. If you were to throw away some of your draws and then take an average, its precision would be no higher than the precision you would expect from a simple random sample from the posterior distribution of size ESS. It is similar to if you had a sample of data that were generated by an AR(1) process and you wanted to estimate the mean of the process, you would use the whole sample rather than some fraction of it.

Ben

[Andrew Gelman]

If you have close to 100% efficiency I doubt you need thousands of draws.

[Michael Betancourt]

PCA is going to be subject to all kinds of noise

and I’ll offer that I don’t think it will be all that

useful.

[Mike Boedigheimer]

@All. Thanks, you have help fill in many gaps for me. Here is my synopsis

For univariate and multivariate parameters. HMC has weak but easily detectable correlations with lags that run into the 10s of simulations (for the above case). These are not revealed by Neff or Rhat. This was not surprising to you but I guess it was to me. 

You can't sample all the tails all the time. In the end, it appears the predictive power of the QDA classifier comes from correlations of threads sampling (by chance) different tail areas. 

The effect is small, and one is probably better off not thinning (at least in this case). As far as I can tell, thinning isn't evil, but just less efficient. You get better estimates by keeping the whole set of simulations than throwing out some to avoid correlations. This doesn't totally sit well with me because I know this leads to highly biased results in repeated measures analysis, which seems somewhat analogous to me. But, I've come to agree it seems to be true here.

I'm not sure how to determine the number of draws I need for any given purpose. I'm very interested in this as well, but that wasn't my question right now. 

Am I on the right page yet? 

@Mike. I agree that PCA isn't the right tool here. However, if you are referring to knowing how to count underlying latent variables, it is actually pretty interesting to me to estimate the distribution of the eigenvalues by permutations. These are compared to the original eigenvalues. 

[Bob Carpenter]

Autocorrelation length depends on model.

Using all draws in a chain is not just a good idea, it's
a theorem that it performs better and is unbiased.

For number of draws you need, look at MCMC standard error on
values you're estimating such as parameters, predictions,
variances, or event probabilities.

The first chapter of the MCMC Handbook is free online and
useful for the basic math of every aspect of this discussion.
The only place we vary is recommending multiple chains to
detect non-convergence issues rather than relying totally on
theoretical asymptotics.

- Bob

[Michael Betancourt]

I’ve said it before I’ll say it again THERE IS NO SUCH THING AS _AN_ 

AUTOCORRELATION IN A MARKOV CHAIN.

Consider a target distribution, pi, a Markov transition, T, and some

function whose expectation we want to take, f.  Asymptotically

\lim_{N \rightarrow \infty} \sum_{n = 1}^{N} f(x_{n}) = E_{pi} [ f ].

Given some nicer regularity conditions (usually we need at least

geometric ergodicity) then the MCMC estimator

\sum_{n = 1}^{N} f(x_{n})

converges in distribution to 

N( E_{pi} [ f ], sigma^{2} / ESS )

where sigma is the true variance of f under pi and ESS is the

effective sample size,

ESS = N / (1 + \sum_{l = 1}^{\infty} rho_{l} ),

where l is the lag-l autocorrelation of f.  This is called the MCMC

central limit theorem.

The important thing to keep in mind here is that the autocorrelations

depend on f!  There is no way to define a canonical autocorrelation

for the chain itself.  So some functions can have large autocorrelations

and some can have small autocorrelations _in the same chain_.  

Intuitively, the more a function varies the more it relies on the fine

details of the target distribution, the longer it will need the Markov

chain to converge, and the larger the autocorrelations will be.

Your PCA is just constructing a function as sensitive to the spatial

scales as possible, yielding larger autocorrelations but that has

absolutely nothing to do with the other expectations and does not

indicate that those ESS calculations are wrong.

The only caveat is that we have access to only ESS _estimates_

and we have to worry about how accurate those estimates can

be.  But these kinds of analyses do not provide any such utility.

Finally, assuming that the MCMC CLT holds you define the number

of steps you need by first determining how well you need to estimate

E_{pi} [ f ] and then running until

MCSE = sqrt{ sigma^{2} / ESS }

achieves that desired precision.  MCSE is printed along with the

other Stan summaries.

[Mike Boedigheimer]

Thanks! I only needed to read the first sentence. 

"A sequence X1, X2, ... of random elements of some set is a Markov chain if the conditional distribution of Xn+1 given X1, ... , Xn depends on Xn only. T"

I'm clearly not dealing with a Markov chain. Because clearly the distribution of Xn+1 depends on Xn and Xn-1. 

[Mike Boedigheimer]

I hope you know I'm kidding. I will need time to digest. 

[Avraham Adler]

Someone has to collect and collate these gems of Michael for display on a website/wiki that supports some kind of LaTeX rendering!

[Krzysztof Sakrejda]

I want to know when the book is coming out.

[Bob Carpenter]

Realistically maybe a year or two. And
we'll likely target beginner's more than Michael's
audience on our mailing list.

Andrew also plans to rewrite the Gelman and
Hill book as two books based entirely on Stan.
That's waiting on someone to code max marginal likelihood,
but right now the algorithm's only pseudocode in
a paper (though we now have the higher-order autodiff
to do it thanks to Rob).

- Bob

[Mike Boedigheimer]

I agree. Once I got passed the shouting, there were all those little gems right there in front of me. I think I've seen them before but there is a gap recognizing the connections between the theory and the problems I'm working on.

So thanks to Michael for "saying it again"

[Michael Betancourt]

I shout only because this has been answered

multiple times on the list.  Please, please, please

consult the archives before asking a high level

question like this.  I’m happy to clarify details or

answer specifics if previous answers aren’t clear,

but at least check those old answers our first.

[Mike Boedigheimer]

I know exactly why you shout and I don't blame you. That is why I was so thankful. I understand that you have to explain these things over and over. It may be hard for you to see that the answers aren't obvious. I looked very hard for posts that discussed observations like those in my OP and couldn't find any. So even though I was pretty sure it had all been answered before, I still didn't get it, and asked the question here.  But, the original responses did not lead me to believe the answers were well know. In the end, it is true that everything I needed was actually contained in other posts, but I was not able to see how (or if) they applied to my situation until you explained it. 

[Michael Betancourt]

Oh, none of this stuff is trivial.  In fact it’s a nightmare

trying to understand the literature from an applied

perspective — treatments are either hugely theoretical

or simplify to the point where they miss important points

with practical consequences.  We’ll hopefully start adding

more and more of this stuff to the manual, but papers

come first!

[Andrew Gelman]

Forward all questions to /dev/ripley/