Stan has RNGs; 1.3.0 release next week?

[Bob Carpenter]

I just exposed Peter's RNG implementationss to the Stan
language and added parsing/compilation tests (Peter already had
tests for the RNGs themsevles). I also added minimal
doc on the changes to the language (RNGs only valid in
generated quantities), but haven't written a new section
on posterior predictive checks (though I did refer readers
to the relevant BDA chapter).

Daniel and I were discussing rolling out a 1.3.0 release
with Daniel's refactoring of math/agrad, these new RNG
implementations, and all of our already-in bug fixes.
Hopefully some time early next week.

Then we can put Michael's refactoring of mcmc, which looks
great, into the next release. (Hope that's OK with everyone.)

Here's the list of available RNGs from the test program; the
commented out ones are still on the to-do list. We still need
to add some domain tests for the RNGs and add them to the doc.
Neither should take too long because they're both straightforward
cut-and-paste jobs (more on refactoring domain tests later).

These RNGs can only be used in the generated quantities
block (at least until we get an initialization block).
Here's what we have:

n <- bernoulli_rng(0.5);
// n <- bernoulli_logit_rng(0.0);
n <- binomial_rng(15,0.3);
// n <- binomial_logit_rng(15,-1.2);
n <- beta_binomial_rng(42, 0.3, 1.9);
// n <- hypergeometric_rng(5,4,9);
n <- neg_binomial_rng(1.2,3.9);
n <- ordered_logistic_rng(1.9,theta);
n <- poisson_rng(2.7);
// n <- poisson_log_rng(2.7);

n <- categorical_rng(theta);
ns <- multinomial_rng(theta,20);

z <- normal_rng(0,1);
// z <- exp_mod_normal_rng(0.2,1.9,1.5);
// z <- skew_normal_rng(1.1, 2.3, 9.8);
z <- student_t_rng(3,1,2);
z <- cauchy_rng(-1.0,2.9);
z <- double_exponential_rng(1.0,3.8);
z <- logistic_rng(1.0,3.8);
// z <- gumbel_rng(-1.0,2.31);
z <- lognormal_rng(-1.0, 3.6);
z <- chi_square_rng(4.1);
z <- inv_chi_square_rng(4.1);
z <- scaled_inv_chi_square_rng(1.0, 3.0);
z <- exponential_rng(2.9);
z <- gamma_rng(1.0, 3.0);
z <- inv_gamma_rng(0.1,0.1);
z <- weibull_rng(1.0,2.0);
z <- pareto_rng(0.1, 1.5);
z <- beta_rng(110.0, 250.1);
z <- uniform_rng(-1.0, 1.0);

theta <- dirichlet_rng(alpha);
v <- multi_normal_rng(mu,Sigma);
// v <- multi_normal_prec_rng(mu,Sigma);
// v <- multi_normal_cholesky_rng(mu,L);
v <- multi_student_t_rng(3.0,mu,Sigma);
Sigma <- wishart_rng(3.0,Sigma);
Sigma <- inv_wishart_rng(3.0,Sigma);
// Sigma <- lkj_corr_rng(2.5);
// Sigma <- lkj_corr_cholesky_rng(3.0);
// DEPRECATE: Sigma <- lkj_cov_rng(Sigma, mu, alpha, 2.0);

- Bob

[Michael Betancourt]

What's the time scale? I should be putting the final touches on gm/command.hpp today, which would make us ready to merge hmc_refactor. Of course, the user impact would be minimal, so I'm happy to push this to the next release and incorporate updated outputs and the like.

[Bob Carpenter]

OK, then!

I thought you still had to implement the command and
all the different samplers and that it wasn't clear
you'd have time before the move. If we can get this in
before you hit the road, that'd be fantastic.

There's been so much change that perhaps we should do a beta
of this release? It'll only be useful if people try it and
give us feedback if it breaks.

- Bob

[Andrew Gelman]

If it's possible, it could be fun to throw in the time-scheduled version of Stan if that's not too hard to do (the version where the user specifies the total running time in seconds or minutes and then Stan divides up the runtime accordingly).
A

On Apr 3, 2013, at 2:19 AM, Bob Carpenter wrote:

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

[Ben Goodrich]

On Wednesday, April 3, 2013 2:19:53 AM UTC-4, Bob Carpenter wrote:

Here's what we have:

   Sigma <- lkj_corr_rng(3,2.5);
   L <- lkj_corr_cholesky_rng(3,3.0);

I added these two. But we still don't have a Choleksy factor primitive type. I hope that can go into the next release.

Unrelatedly, there is this warning

In file included from src/test/math/seq_view_test.cpp:4:
src/stan/math/seq_view.hpp:6:9: warning: 'EIGEN_DENSEBASE_PLUGIN' macro redefined
#define EIGEN_DENSEBASE_PLUGIN "stan/math/EigenDenseBaseAddons.hpp"
        ^
src/stan/math/matrix/Eigen.hpp:4:9: note: previous definition is here
#define EIGEN_DENSEBASE_PLUGIN "stan/math/matrix/EigenDenseBaseAddons.hpp"
        ^

Ben

[Daniel Lee]

Those warnings should be gone in master now.

[Bob Carpenter]

We'll get to that. It's a fair bit of work because
of the need to modify all the I/O, R-hat and n_eff
routines on the Stan and RStan side, which now
expect equal length chains. (Or at least I think they
do. If we can solve that problem, the timing itself
is pretty easy.)

I had some sense from the meeting that Andrew wanted to
weight chains by number of iterations in computing
R-hat and n_eff. We originally weighted each chain
evenly to be more conservative about convergence.
But if we do that, small stuck chains will drag the
whole set of chains down.

- Bob

[Andrew Gelman]

I'm not sure what to do about that. Probably small chains should drag the whole set of chains down. People can of course feel free to use simulations for inference before convergence has been reached.

[Bob Carpenter]

On 4/3/13 5:41 PM, Ben Goodrich wrote:
> On Wednesday, April 3, 2013 2:19:53 AM UTC-4, Bob Carpenter wrote:
>
> Here's what we have:
> Sigma <- lkj_corr_rng(3,2.5);
> L <- lkj_corr_cholesky_rng(3,3.0);
>
>
> I added these two.

Great.

> But we still don't have a Choleksy factor primitive type. I hope that can go into the next release.

Me too.

We're hoping to refactor the code generator and
ast next.

But if you want, we can sit down and do it together.
I'm pretty much out until next Thursday, but
after that, my schedule's wide open.

Treating it the same way as we treat cov_matrix
and corr_matrix will be easy (that is, it just
reads and writes into a matrix). If we want to do
something fancy with the Cholesky factor usage
(like make it into a triangular view where possible),
that's going to be much more work wrestling with both
code generation and the expression templates in Eigen.

> Unrelatedly, there is this warning
>
> |
> In file included from src/test/math/seq_view_test.cpp:4:
> src/stan/math/seq_view.hpp:6:9: warning: 'EIGEN_DENSEBASE_PLUGIN' macro redefined
> #define EIGEN_DENSEBASE_PLUGIN "stan/math/EigenDenseBaseAddons.hpp"
> ^
> src/stan/math/matrix/Eigen.hpp:4:9: note: previous definition is here
> #define EIGEN_DENSEBASE_PLUGIN "stan/math/matrix/EigenDenseBaseAddons.hpp"

Daniel was going to take (or has already taken) care of that one.

- Bob

[Bob Carpenter]

On 4/3/13 6:49 PM, Andrew Gelman wrote:
> I'm not sure what to do about that. Probably small chains should drag the whole set of chains down.

That was the original plan.

Right now, we wait until each chain has run a given
number of iterations. What often happens is that one
chain goes much more slowly than the others, so we have
to wait for it before analyzing results.

I thought part of the idea of timing was that we'd finish
when most of the chains were done and downweight the laggards.
That'd make the whole process much faster to get to a solution
at a given R-hat or n_eff, though it's cheating relative
to the more conservative R-hat and n_eff definitions.

- Bob

[Andrew Gelman]

I think it could work better to run them all for the same length of time?

[Bob Carpenter]

Let me make a concrete plan:

1. User specifies a total run time t.

2. We run warmup until the first iteration ending at a time past t/2.

3. We then run sampling iterations until the first iteration past time t.

4. Analyze ragged R-hat and n_eff with each chain weighted equally.

Is this what you (Andrew) had in mind? If so, that's what
we can implement. I don't think anyone else has a strong
opinion on how this should go.

I'm not sure what this'll do to computation compared to running
each chain the same number of iterations.

Timed setup: If one chain is slow to converge in wall time,
the other chains will keep running for the same amount of time, which
they didn't technically need to do.

# Iterations condition: If one chain takes more iterations to converge,
the other chains will have to keep running for that many iterations,
which they didn't technically need to do.

What we really want is to run until each chain has converged :-)

What I'm seeing in practice is that some chains are slow in wall
time because adaptation gets stuck, but by the time it gets unstuck,
they sample just fine as measured by number of iterations. In
the stochastic volatilit models, it'd often take one chain four or
five times as many iterations to converge as the others (1000 vs.
200, say) in the untransformed model.

- Bob

[Michael Betancourt]

Might be worth looking into sharing adaptation information amongst the chains to help guard against single chains falling into bad adaptations.

[Andrew Gelman]

That makes sense to me. After warm-up, the chains can then be run independently.

[Andrew Gelman]

Actually, let's think more about this. Given that we are optimizing over global parameters, I think we should definitely be sharing info in the adaptation. To start with, this means that when we're running standard serial Stan, we should not be doing all of chain 1, then all of chain 2, etc. One approach would be to thread everything (i.e., 1 iteration of chain 1, 1 iteration of chain 2, 1 iteration of chain 3, etc), another approach would be to break it into blocks (10 of chain 1, 10 of chain 2, etc).

Then how best to share info? The simplest approach would be, when updating each chain (during warm-up) to randomly select the tuning parameters as estimated so far from one of the other chains. Another, more gentle, approach, would be do occasionally swap tuning parameters. We might also want some sort of competition or birth-death rule so that the losers get eliminated.

I'm curious if Matt has any thoughts on this.

[Bob Carpenter]

It's not so much bad adaptations that
are the problem so much as that it sometimes
takes a while to get there, either in time
or in iterations or both.

- Bob

[Andrew Gelman]

Hi--I guess we should think about how we do adaptation before implementing this? It seems like we should switch to a better approach than running the chains separately, one at a time.
A

On Apr 3, 2013, at 7:25 PM, Bob Carpenter wrote:

[Nigel/essence]

In my implementation of NUTS, I do as Andrew suggests (i.e., 1 iteration of chain 1, 1 iteration of chain 2, 1 iteration of chain 3, etc.), and like the idea of selecting adaptation parameters from other chains, but I also throw in some non-NUTS iterations (e.g. DE or RWM), and a resetting of adaptation parameters, to overcome the sticking problem. The sticking problem is definitely an issue. Robust identification of when a chain is sticking may be useful, so corrective action can be taken. The naïve approach would be when the step size collapses, or the number of depths is persistently hitting the limit.

[Nigel/essence]

ps, DE and similar of course require that all chains are run in parallel (i.e., 1 iteration of chain 1, 1 iteration of chain 2, 1 iteration of chain 3, etc.), so this is another side benefit of that architecture.

[Bob Carpenter]

On 4/6/13 10:27 AM, Nigel/essence wrote:
> In my implementation of NUTS, I do as Andrew suggests (i.e., 1 iteration of chain 1, 1 iteration of chain 2, 1 iteration of chain 3, etc.),

Do you try to parallelize, or does this force
everything to run in a single thread/process?

I'm guessing from Andrew's past comments that
he'd want to run multiple instances of a multi-chain
process anyway.

The emcee ("MCMC hammer") Python package implements
Goodman and Weare's ensemble-based sampler by splitting
the walkers in their ensembles into batches and running
the batches in parallel. I think I've said this before, but
the Goodman/Weare method is different than DE/DREAM in the way it
proposes and how it has to compute Metropolis [took me a while
to see it, even at the board with a very patient Jonathan
Goodman].

We'll be adding some of these ensemble methods, including
at least DE and Goodman/Weare's walker method to Stan
right after we get the MCMC base classes refactored a bit.

- Bob

[Nigel/essence]

The obvious thing is to make each chain a separate thread, but in fact in my case (multiple gaussian processes) I use threads within the likelihood calculation, for each gaussian process.