feature/flexible_adaptation tests

[Daniel Lee]

---------- Forwarded message ----------

From: Michael Betancourt <betan...@gmail.com>
Date: Tue, Nov 26, 2013 at 1:22 PM
Subject: feature/flexible_adaptation Tests
To: "stan...@googlegroups.com" <stan...@googlegroups.com>

Acid tests for Volume 1 of the BUGS models.  There are instances where the ESS and ESS/time decrease, but I believe these are due to develop being biased.  For example, notice the behavior of Mice and Rats.  Big reductions of bias in the new branch, as expected.

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

[Daniel Lee]

Thanks for sending along the plots. They're useful, but I think I'm not seeing the improvement branch over develop.

I have a lot of questions regarding the generation of these plots (questions about the branch are lower):

1. How did you generate the plots? I'd like to have Jenkins automatically create something like that if we can.
2. Specifics on how it was generated? Was each model run once? Is each count in the histogram a single parameter? Or are there multiple chains in there?
3. Could you dump out the number of sampling and warmup iterations in the plots? It looks like they might be too low to guarantee convergence.
4. How are you encoding the expected means?
   
5. Are the MCSE taken from the BUGS runs also? Or are they the MCSE of the parameter for the chain?
6. Could we replot as (mean - expected_mean) / expected_mean to get it on an error ratio scale that's comparable?

Questions regarding the branch:

1. The claim is "big reduction of bias in the new branch." I take it you're referring to the first plot on all the pages, the (estimate - mean) / MCSE histogram. If so, I don't see a systematic reduction in bias.
2. Given that it doesn't look like we have dominant behavior, can we independently test the target acceptance rate and new adaptation?
3. Also... if we're only looking at single runs, this might be iffy and maybe we should figure out how to push it to the cluster to do a real comparison of multiple, long runs.
4. Is the new branch really that much slower than develop?

A general comment: in the plots that show both develop and the new branch not recovering the parameters, the R-hats look good, effective sample sizes are high, but the estimates are biased. That's a little worrisome.

I'm going to try to paste in a chart that's my scorecard on the set of runs you've provided. My tally is: 6-16-3-2 (that's flexible_adaptation - develop - draws - losses). This is all hinged on me interpreting the graphs correctly, so if there's something I've missed, please let me know. The last column in the table indicates the branch with problems for that model. Of course, all done by eye.

Model	Less Biased	Faster	Winner	Something's wrong
blocker	flexible_adaptation	develop
bones	draw	develop	develop
dogs	flexible_adaptation	develop
dyes	develop	draw	develop	flexible_adaptation
epil	draw	develop	develop
equiv	flexible_adaptation	develop	flexible_adaptation	develop
inhalers	draw	develop	develop
kidney	flexible_adaptation	flexible_adaptation	flexible_adaptation	develop
leuk	draw	draw	draw	both
leukfr	draw	develop	develop
lsat	develop	draw	draw
magnesium	flexible_adaptation	develop	flexible_adaptation	develop
mice	flexible_adaptation	develop	flexible_adaptation	develop
oxford	develop	develop	develop
pump	draw	draw	draw	both
rats	flexible_adaptation	develop	flexible_adaptation	develop
salm	draw	develop	develop	both
salm2	draw	develop	develop	both
seeds	draw	develop	develop
stacks_a_normal	draw	develop	develop	both
stacks_b_dexp	draw	develop	develop
stacks_c_t4	develop	develop	develop
stacks_d_normal_ridge	develop	develop	develop	flexible_adaptation
stacks_e_dexp_ridge	draw	develop	develop	both
stacks_f_t4_ridge	draw	develop	develop	both
surgical	draw	develop	develop
surgical_simple	flexible_adaptation	develop	flexible_adaptation

[Michael Betancourt]

> • How did you generate the plots? I'd like to have Jenkins automatically create something like that if we can.

Bash scripts, gnuplot, and pdflatex.

> • Specifics on how it was generated? Was each model run once? Is each count in the histogram a single parameter? Or are there multiple chains in there?

Each model is run 50 times, and each histogram is a dump of all parameters for all 50 runs.

> • Could you dump out the number of sampling and warmup iterations in the plots? It looks like they might be too low to guarantee convergence.

Stan was run with defaults, so 1000 warmup and 1000 sampling. This is meant to test the default behavior of Stan (see below).

> • How are you encoding the expected means?

I'm grepping them from the model tests.

> • Are the MCSE taken from the BUGS runs also? Or are they the MCSE of the parameter for the chain?

The MCSE is taken from the Stan runs.

> • Could we replot as (mean - expected_mean) / expected_mean to get it on an error ratio scale that's comparable?

The problem is that that becomes less interpretable. The MCSE standard error really provides the relevant scale for
comparing agreement across models.

> Questions regarding the branch:
>
> • The claim is "big reduction of bias in the new branch." I take it you're referring to the first plot on all the pages, the (estimate - mean) / MCSE histogram. If so, I don't see a systematic reduction in bias.

The point of this branch is to make adaptation far more robust out of the box,
and given that many of the BUGS models are "easy"
there shouldn't be a systematic reduction. The important features to notice are

a) Rhat gets systematically better (indicating faster convergence)

b) In many models outliers vanish (for example, equiv). Note that these
outliers are not evident in the Rhat distributions and so a naive user
would assume the sampler is working fine.

c) In a few models (kidneys, mice, rats) develop looks like it has converged
by the usual metrics but the chains are way off. The new branch, however,
finds the right mode in each case.

> • Given that it doesn't look like we have dominant behavior, can we independently test the target acceptance rate and new adaptation?

We might as well run two branches at that point. Honestly the effects of
the two branches will be similar.

> • Also... if we're only looking at single runs, this might be iffy and maybe we should figure out how to push it to the cluster to do a real comparison of multiple, long runs.

See above.

> • Is the new branch really that much slower than develop?

They can be, but also remember that many of those models run really fast
so the difference in timing isn't really noticeable in practice. Again, one
has to be careful in assuming that develop is getting the right answer
and not sacrificing robustness for speed.

> A general comment: in the plots that show both develop and the new branch not recovering the parameters, the R-hats look good, effective sample sizes are high, but the estimates are biased. That's a little worrisome.

Yes, which is why I'm eager to get these fixes in. Similar behavior is
evident in many hierarchical models which motivated the changes
in the new branch in the first place.

> I'm going to try to paste in a chart that's my scorecard on the set of runs you've provided. My tally is: 6-16-3-2 (that's flexible_adaptation - develop - draws - losses). This is all hinged on me interpreting the graphs correctly, so if there's something I've missed, please let me know. The last column in the table indicates the branch with problems for that model. Of course, all done by eye.

Note that the cases you peg flexible_adaptation with "something's wrong"
are pretty mild. The (mean- expected_mean) / MCSE usually is just
pivoting around 0 but not really getting any worse. Again I have to warn
caution in taking the develop performance too seriously as many of these
models haven't been carefully studied to ensure that they are getting
robust estimates in the first place.

Also recall the use of Stan defaults, not tens of thousands of warmup
interations.

[Bob Carpenter]

Can we change colors from green/red to something more
color-blind friendly? (I hope you don't think I'm
obsessed with color --- I'm much more obsessed with
typography.)

I'd find stacked histograms easier to
interpret than the overlays. I tmay be my slight red-green
color blindness, but I have a hard time trying to "see"
the histograms themselves without a lot of work.

On 11/27/13, 4:24 AM, Michael Betancourt wrote:

>> � Could we replot as (mean - expected_mean) / expected_mean to get it on an error ratio scale that's comparable?

>
> The problem is that that becomes less interpretable. The MCSE standard error really provides the relevant scale for
> comparing agreement across models.

I think they do two different things.

The MCSE plot is more interpretable for a single model in
telling us whether our MCSE estimates are accurate.

I like the (mean - expected)/expected for simply comparing
fits. It puts all the parameters on the same scale, too.

With MCSE, I can't tell which setup actually got closer
to the true values, only which setup thinks it got closer
based on its MCSE estimate.

Can we compromise on having both?

>> Questions regarding the branch:
>>
>> � The claim is "big reduction of bias in the new branch." I take it you're referring to the first plot on all the pages, the (estimate - mean) / MCSE histogram. If so, I don't see a systematic reduction in bias.

>
> The point of this branch is to make adaptation far more robust out of the box,
> and given that many of the BUGS models are "easy"
> there shouldn't be a systematic reduction. The important features to notice are
>
> a) Rhat gets systematically better (indicating faster convergence)

I think the relevant notion of "fast" here is wall time, not R-hat over
1000 iterations.

> b) In many models outliers vanish (for example, equiv). Note that these
> outliers are not evident in the Rhat distributions and so a naive user
> would assume the sampler is working fine.

I'd really really like to see (mean - expected)/expected to see this.
Is the outlier from a bad MCSE estimate or is that estimate really
far away from the true value?

> c) In a few models (kidneys, mice, rats) develop looks like it has converged
> by the usual metrics but the chains are way off. The new branch, however,
> finds the right mode in each case.

Doesn't convergence just mean that it's mixing? And isn't
that consistent with having a high MCSE? (Those
weren't rhetorical questions!)

This is the kind of reason why I want to see (mean - expected)/expected,
to see if we're really converging or not.

>> � Given that it doesn't look like we have dominant behavior, can we independently test the target acceptance rate and new adaptation?

>
> We might as well run two branches at that point. Honestly the effects of
> the two branches will be similar.

I don't understand what Michael means here. What do you mean
by run two branches? I think we can run develop with a target
acceptance of 0.8 and flexible adaptation with an acceptance of 0.5.

Running a bunch of acceptance rates would be even better, to match
the empirical performance with the theory for hierarchical models. But
I know this stuff takes forever to run in bulk.


>> � Is the new branch really that much slower than develop?

>
> They can be, but also remember that many of those models run really fast
> so the difference in timing isn't really noticeable in practice. Again, one
> has to be careful in assuming that develop is getting the right answer
> and not sacrificing robustness for speed.

Agreed. That's why I want (mean - expected)/expected --- to see
if we're getting the right answer!

>> A general comment: in the plots that show both develop and the new branch not recovering the parameters, the R-hats look good, effective sample sizes are high, but the estimates are biased. That's a little worrisome.
>
> Yes, which is why I'm eager to get these fixes in. Similar behavior is
> evident in many hierarchical models which motivated the changes
> in the new branch in the first place.

But you see this in the new branch, too. salm is an example where
flex is biased high and develop biased low, and both have Rhats
under 1.01.

What I'm worried about is that for models like salm, develop
seems to run about twice as fast.

Practically speaking, our users are going to balk if they see
their time per effective sample go out. They're already balking
at 2.0 given that it seems to be slower per iteration.

>> I'm going to try to paste in a chart that's my scorecard on the set of runs you've provided. My tally is: 6-16-3-2 (that's flexible_adaptation - develop - draws - losses). This is all hinged on me interpreting the graphs correctly, so if there's something I've missed, please let me know. The last column in the table indicates the branch with problems for that model. Of course, all done by eye.

I see the chart, but I don't understand the summary.

For speed, it only makes sense to me to compare speed when
both are getting the right answer (again, I want to see
(mean - expected)/expected to judge this). Given the way
we're computing, MCSE and n_eff and R_hat are all closely
related because they both rely on the cross-chain evals.

> Note that the cases you peg flexible_adaptation with "something's wrong"
> are pretty mild. The (mean- expected_mean) / MCSE usually is just
> pivoting around 0 but not really getting any worse. Again I have to warn
> caution in taking the develop performance too seriously as many of these
> models haven't been carefully studied to ensure that they are getting
> robust estimates in the first place.

I think we do know the answer. BUGS ran them for a gazillion
iterations and then we matched that by running many iterations
in Stan, so I think the expected values are pretty reliable.
If they're not reliable, then the (mean - expected)/MCSE is
irrelevant, as is (mean - expected)/expected.

Any idea what's up with the pivoting? It's pretty extreme with
rats and salm and stacks_a_normal, and systematically seems to be
biased high for flex and low for develop.

And just in case you didn't get it earlier, I'd really like to see
(mean - expected)/expected in order for me to understand the scale
of the bias relative to the true values rather than relative to MCSE.

I think the new flex adaptation will actually look better on the
(mean - expected)/expected basis, because it should have lower
MCSE estimates and hence inflated looking bias on the
(mean - expected)/MCSE basis.

dyes and stacks_d_normal_ridge looks much worse in flex for
the (mean - expected)/MCSE basis and leuk looks mildly worse
(it might help to center the histograms around 0 to make this more
visually apparent in the plots). Some look much better in the the
flex branch, kidney and mice in particular.

> Also recall the use of Stan defaults, not tens of thousands of warmup
> interations.

But the Stan defaults are what most people are going to use.

Weren't we going to turn down the number of samples by default?
We have a problem with people just turning Stan on with defaults
and a super-complicated model and then waiting for a weekend for
results.

- Bob

[Michael Betancourt]

> Can we change colors from green/red to something more
> color-blind friendly?

Okay -- just choose two whose sum is pretty (red and blue works well).

> I'd find stacked histograms easier to
> interpret than the overlays. I tmay be my slight red-green
> color blindness, but I have a hard time trying to "see"
> the histograms themselves without a lot of work.

I absolutely hate stacked histograms -- they're almost impossible to read.
I'm happy to tweak the colors but not the histograms.

> I think they do two different things.
>
> The MCSE plot is more interpretable for a single model in
> telling us whether our MCSE estimates are accurate.
>
> I like the (mean - expected)/expected for simply comparing
> fits. It puts all the parameters on the same scale, too.
>
> With MCSE, I can't tell which setup actually got closer
> to the true values, only which setup thinks it got closer
> based on its MCSE estimate.
>
> Can we compromise on having both?

No.

Look, the main issue is that (mean - expected) / expected has no statistical
meaning and makes it impossible to compare between parameters (different
variances, etc).

If you want to discriminate between mean - expected being small vs
MCSE being large then we'd just add a histogram for the MCSE.

>> a) Rhat gets systematically better (indicating faster convergence)
>
> I think the relevant notion of "fast" here is wall time, not R-hat over
> 1000 iterations.

No, it's R-hat/iteration. Develop has a bad habit of getting stuck
in bad tails from which it takes forever to recover (in terms of
number of iterations). flexible_adaptation is specifically designed
to be more robust to this.

>> c) In a few models (kidneys, mice, rats) develop looks like it has converged
>> by the usual metrics but the chains are way off. The new branch, however,
>> finds the right mode in each case.
>
> Doesn't convergence just mean that it's mixing? And isn't
> that consistent with having a high MCSE? (Those
> weren't rhetorical questions!)

Convergence means true mixing, but true mixing isn't always easy to determine.
For example, when latent parameters get stuck (think variances) the rest of
the parameters will still bounce around according to their conditional variances.
It looks like they're sampling fine under the usual diagnostics but they very much
are not.

> This is the kind of reason why I want to see (mean - expected)/expected,
> to see if we're really converging or not.

Converging also means converging in the central limit sense. If the (mean - expected) / MCSE
histograms don't converge to a Gaussian then the chain is boned one way or another.

> I don't understand what Michael means here. What do you mean
> by run two branches? I think we can run develop with a target
> acceptance of 0.8 and flexible adaptation with an acceptance of 0.5.

That's the same as refactoring flexible_adaptation into two branches, which is
more principled then hacking develop.

> Running a bunch of acceptance rates would be even better, to match
> the empirical performance with the theory for hierarchical models. But
> I know this stuff takes forever to run in bulk.

No, these kinds of tests are research topics not Stan tasks. If the hmc_for_hier draft doesn't
convince you that we understand the pathologies then I'm not sure what will.

> Agreed. That's why I want (mean - expected)/expected --- to see
> if we're getting the right answer!

Again, (mean - expected) / expected doesn't tell you anything rigorous.
Only (mean - expected) / MCSE has any proper relevance to MCMC.

> But you see this in the new branch, too. salm is an example where
> flex is biased high and develop biased low, and both have Rhats
> under 1.01.

Those are small biases -- one isn't any better than the other. I'm
referring to the HUGE biases that are fixed under the new branch.

> What I'm worried about is that for models like salm, develop
> seems to run about twice as fast.
>
> Practically speaking, our users are going to balk if they see
> their time per effective sample go out. They're already balking
> at 2.0 given that it seems to be slower per iteration.

No one is doing any rigorous tests to make any proper statement
so I'm not particularly concerned. flexible_adaptation is slower
because the tree_depths are larger. There is some possibility
that the variances are being underestimated, but it's just as likely
that flexible_adaptation is exploring more (which is also more
likely given the larger ESS for that model).

> I think we do know the answer. BUGS ran them for a gazillion
> iterations and then we matched that by running many iterations
> in Stan, so I think the expected values are pretty reliable.
> If they're not reliable, then the (mean - expected)/MCSE is
> irrelevant, as is (mean - expected)/expected.
>
> Any idea what's up with the pivoting? It's pretty extreme with
> rats and salm and stacks_a_normal, and systematically seems to be
> biased high for flex and low for develop.

Running BUGS for a gazillion iterations does not guarantee the
right answer. The hierarchical models that are causing problems
having exactly the geometry that invalidates the usual ergodicity
proofs for RWM and Gibbs.

One possibility with the pivoting (rats isn't a pivot, by the way,
it's just the new branch getting the right answer; I'm referring
to pivoting as the new branch overshooting 0) is that the baseline
values are biased (they couldn't reach far enough into a tail,
for example) and the new branch is going all the way. This
kind of behavior is easy to see for the funnel.

> I think the new flex adaptation will actually look better on the
> (mean - expected)/expected basis, because it should have lower
> MCSE estimates and hence inflated looking bias on the
> (mean - expected)/MCSE basis.

Huh? MCSE only scales the distribution -- it won't correct for any
biases.

>
> dyes and stacks_d_normal_ridge looks much worse in flex for
> the (mean - expected)/MCSE basis and leuk looks mildly worse
> (it might help to center the histograms around 0 to make this more
> visually apparent in the plots). Some look much better in the the
> flex branch, kidney and mice in particular.

Yes, dyes and stacks_d_normal_ridge look suspicious and warrant
further investigation.

> But the Stan defaults are what most people are going to use.

Which is why robust adaptation should focus on performance
given those defaults.

> Weren't we going to turn down the number of samples by default?
> We have a problem with people just turning Stan on with defaults
> and a super-complicated model and then waiting for a weekend for
> results.

I think the timing warning will help there.

[Daniel Lee]

On Wed, Nov 27, 2013 at 11:47 AM, Bob Carpenter <ca...@alias-i.com> wrote:

Can we change colors from green/red to something more
color-blind friendly?  (I hope you don't think I'm
obsessed with color --- I'm much more obsessed with
typography.)

I'd find stacked histograms easier to
interpret than the overlays.  I tmay be my slight red-green
color blindness, but I have a hard time trying to "see"
the histograms themselves without a lot of work.

On 11/27/13, 4:24 AM, Michael Betancourt wrote:

        • Could we replot as (mean - expected_mean) / expected_mean to get it on an error ratio scale that's comparable?

The problem is that that becomes less interpretable.  The MCSE standard error really provides the relevant scale for
comparing agreement across models.

I think they do two different things.

The MCSE plot is more interpretable for a single model in
telling us whether our MCSE estimates are accurate.

I like the (mean - expected)/expected for simply comparing
fits.  It puts all the parameters on the same scale, too.

With MCSE, I can't tell which setup actually got closer
to the true values, only which setup thinks it got closer
based on its MCSE estimate.

Can we compromise on having both?

Questions regarding the branch:

        • The claim is "big reduction of bias in the new branch." I take it you're referring to the first plot on all the pages, the (estimate - mean) / MCSE histogram. If so, I don't see a systematic reduction in bias.

The point of this branch is to make adaptation far more robust out of the box,
and given that many of the BUGS models are "easy"
there shouldn't be a systematic reduction.  The important features to notice are

a) Rhat gets systematically better (indicating faster convergence)

I think the relevant notion of "fast" here is wall time, not R-hat over
1000 iterations.

b) In many models outliers vanish (for example, equiv).  Note that these
outliers are not evident in the Rhat distributions and so a naive user
would assume the sampler is working fine.

I'd really really like to see (mean - expected)/expected to see this.
Is the outlier from a bad MCSE estimate or is that estimate really
far away from the true value?

c) In a few models (kidneys, mice, rats) develop looks like it has converged
by the usual metrics but the chains are way off.  The new branch, however,
finds the right mode in each case.

Doesn't convergence just mean that it's mixing?  And isn't
that consistent with having a high MCSE?  (Those
weren't rhetorical questions!)

This is the kind of reason why I want to see (mean - expected)/expected,
to see if we're really converging or not.

        • Given that it doesn't look like we have dominant behavior, can we independently test the target acceptance rate and new adaptation?

We might as well run two branches at that point.  Honestly the effects of
the two branches will be similar.

I don't understand what Michael means here.  What do you mean
by run two branches?  I think we can run develop with a target
acceptance of 0.8 and flexible adaptation with an acceptance of 0.5.

Running a bunch of acceptance rates would be even better, to match
the empirical performance with the theory for hierarchical models.  But
I know this stuff takes forever to run in bulk.

        • Is the new branch really that much slower than develop?

I was just eyeballing things. If one clearly had less bias, it won. If that looked about the same, the faster one won. If both were even for bias and speed, it was a draw. There were a few cases where both clearly lost against the example, and then I didn't make a judgement.

I think this is worth repeating. Although the new branch fixes some of these issues, there are some models for which the new branch is broken. 'stacks_d_normal_ridge' is an example of where the new branch is sampling in the wrong place, not just equally biased.

 

Also recall the use of Stan defaults, not tens of thousands of warmup
interations.

But the Stan defaults are what most people are going to use.

Weren't we going to turn down the number of samples by default?
We have a problem with people just turning Stan on with defaults
and a super-complicated model and then waiting for a weekend for
results.

- 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.

[Daniel Lee]

On Wed, Nov 27, 2013 at 1:13 PM, Michael Betancourt <betan...@gmail.com> wrote:

> Can we change colors from green/red to something more
> color-blind friendly?

Okay -- just choose two whose sum is pretty (red and blue works well).

I'll stay far away from this discussion.

 

> I'd find stacked histograms easier to
> interpret than the overlays.  I tmay be my slight red-green
> color blindness, but I have a hard time trying to "see"
> the histograms themselves without a lot of work.

I absolutely hate stacked histograms -- they're almost impossible to read.
I'm happy to tweak the colors but not the histograms.

The overlays, in this case, are hard to read. They might have to do with color choice, so if you fix it, hopefully it'll be clearer.

I don't have any trouble with stacked histograms, as long as the axes are on the same scale.

 

> I think they do two different things.
>
> The MCSE plot is more interpretable for a single model in
> telling us whether our MCSE estimates are accurate.
>
> I like the (mean - expected)/expected for simply comparing
> fits.  It puts all the parameters on the same scale, too.
>
> With MCSE, I can't tell which setup actually got closer
> to the true values, only which setup thinks it got closer
> based on its MCSE estimate.
>
> Can we compromise on having both?

No.

Look, the main issue is that (mean - expected) / expected has no statistical
meaning and makes it impossible to compare between parameters (different
variances, etc).

If you want to discriminate between mean - expected being small vs
MCSE being large then we'd just add a histogram for the MCSE.

Can you clarify the use of MCSE? I'm assuming the MCSE used was the estimate for each parameter within each run. So each is normalized by its own MCSE. Not the reference MCSE that's from BUGS or some other consistent normalizing constant between runs and between branches.

One issue that we're going to see is for a particular parameter within a particular run, if the R-hat looks good, we're going to have a high effective sample size, even if the parameter estimate is biased. This is going to shrink things on that graph artificially, where as I think both Bob and I would have a better intuition of the sampler's errors in these parameter estimates if it was relative error.

 

>> a) Rhat gets systematically better (indicating faster convergence)
>
> I think the relevant notion of "fast" here is wall time, not R-hat over
> 1000 iterations.

No, it's R-hat/iteration.  Develop has a bad habit of getting stuck
in bad tails from which it takes forever to recover (in terms of
number of iterations).  flexible_adaptation is specifically designed
to be more robust to this.

Is magnesium one of those cases? I'd be happy to see a particular example that we can work through and really see the differences between develop and the new stuff.

Just to be clear, you're just talking about stuck chains, right? How far off are the parameter estimates? Is this a problem because the n_eff is really low (like 0), and the MCSE is really low? Or are the parameter values just way off?

 

--
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]

[Michael Betancourt]

> Can you clarify the use of MCSE? I'm assuming the MCSE used was the estimate for each parameter within each run. So each is normalized by its own MCSE. Not the reference MCSE that's from BUGS or some other consistent normalizing constant between runs and between branches.

Yes, I'm using the MCSE from each run.

Here's the deal -- the whole utility of MCMC boils down to the Monte Carlo Central Limit Theorem:
assuming a well-behaved sampler and target distribution then the statistic

(mean - expected) / MCSE

converges to a N(0, 1) distribution.

Plotting something like (mean - expected) / expected doesn't mean anything because the scale
has no particularly meaning. The distribution of that statistic will change depending on the number
of samples, the autocorrelation, etc, and will not in general look gaussian. In particular, how do you
determine whether the bias is significant or not?

By plotting (mean - expected) / MCSE for each variable we should be adding a bunch of N(0, 1)
gaussians together, which is another N(0, 1) gaussian, making it easy to interpret the histogram
and any deviations from the expected behavior.

Note that the MC CLT does not always hold for the models we're considered! In particular be wary
of Random Walk Metropolis and the Gibbs sampler for hierarchical models with centered
parameterizations or heavy-tailed distributions. This is one reason why I'm not wedded to
the current "expected" values.

> One issue that we're going to see is for a particular parameter within a particular run, if the R-hat looks good, we're going to have a high effective sample size, even if the parameter estimate is biased. This is going to shrink things on that graph artificially, where as I think both Bob and I would have a better intuition of the sampler's errors in these parameter estimates if it was relative error.

But that's exactly the desired behavior. You want a large ESS to amplify the bias.
I'm not sure what you expect to see without that amplification -- even if the relative
error is small there's something serious going on.

> No, it's R-hat/iteration. Develop has a bad habit of getting stuck
> in bad tails from which it takes forever to recover (in terms of
> number of iterations). flexible_adaptation is specifically designed
> to be more robust to this.
>
> Is magnesium one of those cases? I'd be happy to see a particular example that we can work through and really see the differences between develop and the new stuff.

Yeah, Magnesium is exhibiting that behavior although the statistics are small.

> Just to be clear, you're just talking about stuck chains, right? How far off are the parameter estimates? Is this a problem because the n_eff is really low (like 0), and the MCSE is really low? Or are the parameter values just way off?

Actually it's the other way around. The chain gets stuck at a bad value (say because of overaggressive adaptation)
and wiggles slightly but not according to its true distribution (probably a conditional one). The ESS calculation sees enough variation
to return a large value, and hence small MCSE, even though the samples are biased.

[Bob Carpenter]

On 11/27/13, 4:20 PM, Michael Betancourt wrote:
>> Can you clarify the use of MCSE? I'm assuming the MCSE used was the estimate for each parameter within each run. So each is normalized by its own MCSE. Not the reference MCSE that's from BUGS or some other consistent normalizing constant between runs and between branches.
>
> Yes, I'm using the MCSE from each run.

And for each parameter, because we get different
MCSE's per parameter based on our estimates of posterior
sd and ESS (aka n_eff).

> Here's the deal -- the whole utility of MCMC boils down to the Monte Carlo Central Limit Theorem:
> assuming a well-behaved sampler and target distribution then the statistic
>
> (mean - expected) / MCSE
>
> converges to a N(0, 1) distribution.

Right. It's absolutely what you need to understand the
uncertainty in an estimate based on the posterior mean.

> Plotting something like (mean - expected) / expected doesn't mean anything because the scale
> has no particularly meaning.

What I really intended was

(mean - expected) / abs(expected)

and the meaning is error relative to the true value. With the
abs on the outside (i.e., ignoring bias), it's what Wikipedia
calls "approximation error" and defines here:

http://en.wikipedia.org/wiki/Approximation_error#Definitions

I'd also be happy looking at absolute relative error without
the sign to indicate bias:

abs((estimate - true) / true)

I'm not sure why you say that errors aren't statistical. My
understanding is that abs(true - expected) is absolute error
and (true - expected)^2 is squared error and that these
concepts are widely used in stats. Medians minimize the expectation
of the former and means the latter, for example, and Berger discusses
them at length in his Bayesian decision theory book.

Dividing by abs(expected) just normalizes the scales.

I realize I may be fundamentally confused about something, but I'm not
sure what I'm confused about, since Michael's answers seem to be about
different questions (which I may also not understand --- one problem
of being a novice in a field is not recognizing answers as answers).

> The distribution of that statistic will change depending on the number
> of samples, the autocorrelation, etc, and will not in general look gaussian. In particular, how do you
> determine whether the bias is significant or not?

I'm not trying to determine that. I just want to know which approach
got closer to the expected values (which themselves may be bogus as
Michael notes below).

> By plotting (mean - expected) / MCSE for each variable we should be adding a bunch of N(0, 1)
> gaussians together, which is another N(0, 1) gaussian, making it easy to interpret the histogram
> and any deviations from the expected behavior.

I see why this is useful. But when we're comparing two samplers,
it brings in the fact that we're estimating posterior sd and
posterior ESS (aka n_eff) in both cases. If I understand other points Michael
is making, these calculations are suspect for the same reason as R-hat.

The relative error is not subject to bogus n_eff and posterior sd calcs,
but I do see how it is subject to bogus expected values.

> Note that the MC CLT does not always hold for the models we're considered! In particular be wary
> of Random Walk Metropolis and the Gibbs sampler for hierarchical models with centered
> parameterizations or heavy-tailed distributions. This is one reason why I'm not wedded to
> the current "expected" values.

Understood.

My intention isn't to turn this into a huge effort, but it seems
to call out for using simulated values where we KNOW the true value.

Maybe in that (n-th + 1) paper Andrew is suggesting.

>> One issue that we're going to see is for a particular parameter within a particular run, if the R-hat looks good, we're going to have a high effective sample size, even if the parameter estimate is biased. This is going to shrink things on that graph artificially, where as I think both Bob and I would have a better intuition of the sampler's errors in these parameter estimates if it was relative error.
>
> But that's exactly the desired behavior. You want a large ESS to amplify the bias.
> I'm not sure what you expect to see without that amplification -- even if the relative
> error is small there's something serious going on.

>> No, it's R-hat/iteration. Develop has a bad habit of getting stuck
>> in bad tails from which it takes forever to recover (in terms of
>> number of iterations). flexible_adaptation is specifically designed
>> to be more robust to this.

R-hat/iteration ignores the actual time taken for the iterations. So it's
good in theory to compare algorithms independently of their implementations,
but as a user, I care about time/ESS of the actual implementations.

>> Is magnesium one of those cases? I'd be happy to see a particular example that we can work through and really see the differences between develop and the new stuff.
>
> Yeah, Magnesium is exhibiting that behavior although the statistics are small.
>
>> Just to be clear, you're just talking about stuck chains, right? How far off are the parameter estimates? Is this a problem because the n_eff is really low (like 0), and the MCSE is really low? Or are the parameter values just way off?

That's why I'd like to see relative error.

> Actually it's the other way around. The chain gets stuck at a bad value (say because of overaggressive adaptation)
> and wiggles slightly but not according to its true distribution (probably a conditional one). The ESS calculation sees enough variation
> to return a large value, and hence small MCSE, even though the samples are biased.

I agree that this is the really dangerous case --- when our diagnostics think
we've converged and have a high ESS, but in fact we haven't.

- Bob

[Michael Betancourt]

The problem with relative and absolute errors is a lack of calibration.

Say we run and have (mean - expected) = 1, or (mean - expected) / abs(expected) = 0.1.
Is that small? Is that big? There's no justification in any case so any claim of
"small" or "big" errors has no motivation. The relative error _is_ scale invariant
(although the mean is not invariant to nonlinear transformations so scale invariance
isn't all that helpful anyways) but there's no calibration. Parenthetically,
decision theoretical arguments aren't super helpful here because they
also don't provide helpful scales; this is why isn't so hard to develop a rigorous
criticism of a model in isolation.

On the other hand, since (mean - expected) / MCSE is N(0, 1) so we do have a calibration
and know what is small and what is large.

Now, what pathologies does MCSE introduce? In the worse case ESS = N and we
get a MCSE that's too large and an amplification of any bias. That's a pretty awesome
worse case for diagnostics.

Regarding the validity of the expected values, we have to be careful. This is one
of the reasons why I keep yelling about not taking performance tests too seriously
right now. We don't know if our baseline is correct, and we shouldn't have too
much faith if the baseline was computing in BUGS/JAGS. So please excuse my
frustration when getting a new feature in requires going through this whole mess.

[Daniel Lee]

On Wed, Nov 27, 2013 at 5:27 PM, Michael Betancourt <betan...@gmail.com> wrote:

The problem with relative and absolute errors is a lack of calibration.

Say we run and have (mean - expected) = 1, or (mean - expected) / abs(expected) = 0.1.
Is that small?  Is that big? There's no justification in any case so any claim of
"small" or "big" errors has no motivation.  The relative error _is_ scale invariant
(although the mean is not invariant to nonlinear transformations so scale invariance
isn't all that helpful anyways) but there's no calibration.  Parenthetically,
decision theoretical arguments aren't super helpful here because they
also don't provide helpful scales; this is why isn't so hard to develop a rigorous
criticism of a model in isolation.

On the other hand, since (mean - expected) / MCSE is N(0, 1) so we do have a calibration
and know what is small and what is large.

Now, what pathologies does MCSE introduce? In the worse case ESS = N and we
get a MCSE that's too large and an amplification of any bias.  That's a pretty awesome
worse case for diagnostics.

I think the worst case is the variance of the samples is 0 or close to 0, which happens with chains that are stuck. That completely blows the scale of the (mean - expected) / MCSE out of anything reasonable. I'm not sure what our R-hat calculations say when all samples are the same. I think it might return 1.

 

Regarding the validity of the expected values, we have to be careful.  This is one
of the reasons why I keep yelling about not taking performance tests too seriously
right now.  We don't know if our baseline is correct, and we shouldn't have too
much faith if the baseline was computing in BUGS/JAGS.  So please excuse my
frustration when getting a new feature in requires going through this whole mess.

[Michael Betancourt]

Reran the acid tests and introduced more info into the output.

Note that the variances are the same in each case and the step sizes consistent with the increased target acceptance probability, so there's nothing crazy going on with the actual learning of parameters.  The MCSE also changes as expected from the estimated variance being the same and the ESS increasing.

Spot checks on the problem models are inconsistent with the visual output.  I am investigating as I can.

[Michael Betancourt]

Indexing error -- this is more like it.

[Bob Carpenter]

Thanks for the new colors --- much easier for me to pick out
the individual histograms.

Are the MCSE plots absolute or normalized by expected parameter
size? I ask because I'm surprised there isn't more range given the wide
range of param values. Or are all these models standardized in some way?

What is the plot labeled "Variances"?

If we don't trust the expected values, how am I supposed to
interpret (mean - expected)/MCSE? I see no reason why the
flex adaptation wouldn't be doing better because of the lower
step size. Can I trust the MCSE values? Good cases for
comparison are oxford and salm/salm2, where MCSE looks better for the flex
adaptation but (mean - expected)/MCSE looks better for develop.

On to models. Here's my scorecard for this bout:

SPEED ACCURACY
blocker: dev a bit faster tie
bones: dev 50% faster tie
dogs: dev 50% faster tie
dyes: dev 50% faster tie
epil: dev 100% faster tie
equiv: dev a bit faster tie
inhalers: tie tie
kidney: tie tie
leuk: dev a bit faster tie
leukfr: tie dev a couple bad outliers
lsat: dev a bit slower dev a bit less accurate
magnesium: dev a bit faster tie
mice: dev 50% faster tie
oxford: dev a bit slower dev a bit more accurate
pump: dev 25% faster dev a tiny bit more accurate
rats: dev 100% faster dev a bit more accurate
salm: dev 50% faster dev a bit less accurate
salm2: dev 50% faster dev a bit less accurate
seeds: dev 50% faster dev a bit less accurate
stacks_a: I think the plots are messed up --- there are four colors
stacks_b: dev a bit faster dev a bit more accurate
stacks_c: dev a bit faster dev a bit more accurate
stacks_d: dev a bit faster dev a bit less accurate
stacks_e: dev a bit faster dev a bit more accurate
stacks_f: dev a bit faster dev mostly a bit more accurate except on
a couple instances
surgical: dev 50% faster tied
surg_sim: dev 50% faster dev a bit less accurate


What am I still missing here? Overall, this looks like a clear
win for the develop branch.

It seems that the speed issues are directly related to the
lower step sizes and increased tree depth. This is why I'd
like to see flexible adaptation tested with the old target acceptance
rate. Or the develop branch tested with a higher target acceptance
rate. Can't that be done by just changing the config for acceptance
rates on develop and flexible adaptation branches?

Or maybe we need to test on some other models? As Michael said, these
models are all relatively easy to fit.

- Bob

On 11/28/13, 12:53 PM, Michael Betancourt wrote:
> Indexing error -- this is more like it.
>
>
>
>
>
> On Thu, Nov 28, 2013 at 5:10 PM, Michael Betancourt <betan...@gmail.com <mailto:betan...@gmail.com>> wrote:
>
> Reran the acid tests and introduced more info into the output.
>
> Note that the variances are the same in each case and the step sizes consistent with the increased target acceptance
> probability, so there's nothing crazy going on with the actual learning of parameters. The MCSE also changes as
> expected from the estimated variance being the same and the ESS increasing.
>
> Spot checks on the problem models are inconsistent with the visual output. I am investigating as I can.
>
>

> On Wed, Nov 27, 2013 at 10:36 PM, Daniel Lee <bea...@alum.mit.edu <mailto:bea...@alum.mit.edu>> wrote:

> stan-dev+u...@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

> For more options, visit https://groups.google.com/groups/opt_out.
>
>
> --
> 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 <mailto:stan-dev%2Bunsu...@googlegroups.com>.

[Michael Betancourt]

> Are the MCSE plots absolute or normalized by expected parameter
> size? I ask because I'm surprised there isn't more range given the wide
> range of param values. Or are all these models standardized in some way?

Absolute. The parameters plotting are those with expected values
in the tests so any standardization is incidental.

> What is the plot labeled "Variances"?

Those are the marginal variances used for the mass matrix elements.
The fact that they're the same indicates that the adaptation parameters
are not changing.

> If we don't trust the expected values, how am I supposed to
> interpret (mean - expected)/MCSE? I see no reason why the
> flex adaptation wouldn't be doing better because of the lower
> step size. Can I trust the MCSE values?

Yes, there are no evident pathologies to suggest otherwise.

> Good cases for
> comparison are oxford and salm/salm2, where MCSE looks better for the flex
> adaptation but (mean - expected)/MCSE looks better for develop.

There's really no statistical difference in those distributions, they're
both equally consistent with N(0, 1). In fast, that's true for just about
all of these. Trying to make any judgement by eye is a terrible idea,
these plots are just about identifying gross features.

> It seems that the speed issues are directly related to the
> lower step sizes and increased tree depth. This is why I'd
> like to see flexible adaptation tested with the old target acceptance
> rate. Or the develop branch tested with a higher target acceptance
> rate. Can't that be done by just changing the config for acceptance
> rates on develop and flexible adaptation branches?

The speed different is really due to subtleties in NUTS. As the step size
drops the ESS gets better, but not enough to compensate for the super linear
cost of building the NUTS tree.

> Or maybe we need to test on some other models? As Michael said, these
> models are all relatively easy to fit.

These tests were only to validate that the new code wasn't messing anything
up, not to motivate its existence. If every algorithmic change gets micromanaged
like this we're going to very quickly run end up with stagnant samplers.