almost ready for V1 release

[Bob Carpenter]

The only two things left on the must-do list are:

1. NUTS II
- estimate diagonal mass matrix during warmup

2. Install documentation
- Windows, Linux, and Mac

We should talk about priorities for V2 in the
meeting Wednesday. The only hard-and-fast
top-level priorities are:

3. Sample analysis
- ESS, split R-hat, smallest posterior interval (SPiN),
quantiles, means

4. R interface
- Rcpp

Other things we've talked about are in the top-level
TO-DO.txt file.

- Bob

[Marcus Brubaker]

Hey,

Great news that Stan is so close to a 1.0 release!  Sorry I've been quiet of late, I was travelling (AISTATS, talking up Stan when I could!) and then attempting to catch up from travelling.  Anyway. there were two things on my mental list to look into that would be nice to address for 1.0.

First, the way chains are initialized right now is pretty inflexible.  Basically, the random initialization is uniform on [-1,1] (or maybe [0,1]) which isn't ideal for some models and may not be well dispersed for Rhat statistiics with multiple chains.  Currently the only way around this I've found is to construct initialization files by hand.  For 1.0 it would be nice if the initialization was a little better dispersed and maybe a bit more flexible.  E.G., make the random initialization Normal(0,sigma) where sigma could be set from the command line.  Even better would be Student-t(0,sigma,d) but that's minor.  Ideally, maybe for 2.0, there would be an optional way to specify in the model file an initialization section which could do data and model specific initialization.

The other thing was an easy way to restart chains.  Last I checked this wasn't supported but seems like a pretty important thing to be able to do.  Even just adding an ability to output a text file with the last sample which was in the format of an initialization file would be a pretty good improvement.

Anyway, those are relatively small in the grand scheme of things but would definitely be noticed by people attempting to use it.

Cheers,
Marcus

[Bob Carpenter]

On 5/8/12 7:09 PM, Marcus Brubaker wrote:

> First, the way chains are initialized right now is pretty inflexible. Basically, the random initialization is uniform
> on [-1,1] (or maybe [0,1]) which isn't ideal for some models and may not be well dispersed for Rhat statistiics with
> multiple chains.

Right now, it's uniform on (-2,2).
(line 311 of gm/command.hpp now that
I looked it up).

> Currently the only way around this I've found is to construct initialization files by hand.

That's right. The initial values are
on the model's scale, not the transformed
unconstrained scale.

> For 1.0
> it would be nice if the initialization was a little better dispersed and maybe a bit more flexible. E.G., make the
> random initialization Normal(0,sigma) where sigma could be set from the command line. Even better would be
> Student-t(0,sigma,d) but that's minor.

Either a normal or Student-t would be just as easy
as the uniform.

What I'm worried about is that with very diffuse priors,
the initial values are going to be too far out. For instance,
someone using a uniform(0,1000) parameter on a scale
parameter would get initial values around 500!

> Ideally, maybe for 2.0, there would be an optional way to specify in the model
> file an initialization section which could do data and model specific initialization.

I like that idea.

It should be pretty easy to code with just another
block I would run like the parameters block already.
I'll have to think a bit about constraints and error
checking. Presumably it would also be on the model's
scale, not the unconstrained scale.

What we will need is all the random variate generators.
We also need those for the generated quantities block,
because Andrew wants to use it for "forward sampling"
in a model checking context (it's forward in the
directed graphical model sense of taking fixed values
of parameters from the root down to some point and
then randomly generating the rest -- you don't need
any fancy computation, just generators from the basic
distributions).

One issue is what the syntax should be. Andrew was
thinking of overloading the sampling ~ to work as
a random variate generator. So if you wrote

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

in the derived quantities block (or the parameter
initialization block), you'd get something equivalent
to

beta[n] <- random_normal(0,10);

> The other thing was an easy way to restart chains. Last I checked this wasn't supported but seems like a pretty
> important thing to be able to do. Even just adding an ability to output a text file with the last sample which was in
> the format of an initialization file would be a pretty good improvement.

That's already on the to-do list. It's still not
supported. I struggled with the I/O format. I figured
the dump format for input is too structured for most
people's output needs. But it's not hard to dump our
variables in that format, though so far there's only
a reader, not a writer for the dump format.

It could also stand to be made a bit more robust to
the crazy behavior of R (writing 2:3 instead of c(2,3)
if it can for dimensions, for instance).

> Anyway, those are relatively small in the grand scheme of things but would definitely be noticed by people attempting to
> use it.

Indeed. Thanks for the suggestions.

- Bob

[Marcus Brubaker]

On Tue, May 8, 2012 at 7:54 PM, Bob Carpenter <ca...@alias-i.com> wrote:

> For 1.0

it would be nice if the initialization was a little better dispersed and maybe a bit more flexible.  E.G., make the
random initialization Normal(0,sigma) where sigma could be set from the command line.  Even better would be
Student-t(0,sigma,d) but that's minor.

Either a normal or Student-t would be just as easy
as the uniform.

What I'm worried about is that with very diffuse priors,
the initial values are going to be too far out.  For instance,
someone using a uniform(0,1000) parameter on a scale
parameter would get initial values around 500!

It's true, overly diffuse starts can cause numerical issues.  Particularly with the way we're doing differentiation it's possible for overly diffused initializations to put samples in places with little or no gradient information which would cause very slow sampling.  However, that's always going to be true to one extent or another.

At the end of the day what constitutes a reasonable but diffuse random initialization is always going to be model (if not data) specific. Anything completely generic will be necessarily wrong most of the time :)
 

I was thinking of exactly the same syntax, as it's pretty clear how it works and is consistent with the rest of the language.  Another issue is how to handle unspecified variables and forward references.  E.G., in

paramaters {
  real a;
  real b;
  real c;
}

initialization {
  a ~ normal(b,10)
}

How is c initialized and what is the mean of a?  You could try to be clever and check that people don't do forward references, but that seems needlessly complicated.  One reasonable way is to initialize things using the default method (uniform, normal, whatever) before calling the code implied by this method.

Anyway, when it comes time to do this I have C implementations of random number generators for several distributions (e.g., normal, truncated normal, gamma, beta, beta prime) which I would be happy to contribute to the cause.

The other thing was an easy way to restart chains.  Last I checked this wasn't supported but seems like a pretty
important thing to be able to do.  Even just adding an ability to output a text file with the last sample which was in
the format of an initialization file would be a pretty good improvement.

That's already on the to-do list.  It's still not
supported.  I struggled with the I/O format.  I figured
the dump format for input is too structured for most
people's output needs.  But it's not hard to dump our
variables in that format, though so far there's only
a reader, not a writer for the dump format.

I actually think a raw binary dump of the parameter vector would be a workable solution here, as it shouldn't be used for anything other than restarting the chain.  If you wanted to be a little more portable and fancy, you could use the boost serialization which supports both binary and text formats for the file.  It would also be nice if this file was written after every step of the chain in case execution is terminated abruptly.
 

It could also stand to be made a bit more robust to
the crazy behavior of R (writing 2:3 instead of c(2,3)
if it can for dimensions, for instance).

I have to admit, the R format kinda drives me nuts and right now there are a few things which are a bit problematic.  E.G.,

data {
  real a[5,6];
  matrix(5,6) b;
}

require different data orderings so that a(i,j) == b(i,j) if I recall correctly.  a is row-wise and b is column-wise if memory serves.  Probably something for post-1.0 though since most people seem happy with the R format.
 
Cheers,
Marcus

[Andrew Gelman]

On May 8, 2012, at 8:25 PM, Marcus Brubaker wrote:

> At the end of the day what constitutes a reasonable but diffuse random initialization is always going to be model (if not data) specific. Anything completely generic will be necessarily wrong most of the time :)

Yes, but it's a good general principle to parameterize models so that parameters are of order 1.

>
> initialization {
> a ~ normal(b,10)
> }

I think people will like this sort of thing. Defaults are great, but as a user I like to have control over the settings. One of the things that frustrated me about Bugs/Jags was that I couldn't tell it how to do the sampling.

>
> I actually think a raw binary dump of the parameter vector would be a workable solution here, as it shouldn't be used for anything other than restarting the chain.

I recommend avoiding binaries! I'd prefer to have the object that is R-readable.

[Bob Carpenter]

Andrew -- you will have full access to the
values in R. This is just to make restarting
convenient and efficient.

A binary dump to restart of the unconstrained
parameters would definitely be easiest -- assuming binary
dumps in C++ are at all easy. We could also
just write it out and scan it back in. It
should hardly be the bottleneck to anything.

But it's not going to be much use in R without
the inverse transform back to the model parameter
space.

- Bob

[Marcus Brubaker]

On Tue, May 8, 2012 at 8:52 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

On May 8, 2012, at 8:25 PM, Marcus Brubaker wrote:

> At the end of the day what constitutes a reasonable but diffuse random initialization is always going to be model (if not data) specific. Anything completely generic will be necessarily wrong most of the time :)

Yes, but it's a good general principle to parameterize models so that parameters are of order 1.

True, but sometimes that's hard to do unless you know a fair bit about the data in the first place.
 

>
> initialization {
>   a ~ normal(b,10)
> }

I think people will like this sort of thing.  Defaults are great, but as a user I like to have control over the settings.  One of the things that frustrated me about Bugs/Jags was that I couldn't tell it how to do the sampling.

Agreed, I think this is the ideal solution.  Good, sensible defaults that work well 90% of the time are great but if you can't control things for that other 10% of the time, it can be really frustrating.
 

>
> I actually think a raw binary dump of the parameter vector would be a workable solution here, as it shouldn't be used for anything other than restarting the chain.

I recommend avoiding binaries!  I'd prefer to have the object that is R-readable.

As Bob noted, this would be only for restarting a chain that had finished or been interrupted.  A binary blob isn't ideal but it is the easiest and fastest to implement.  Now that I think about it, a more structured data format may be needed though, since it would be nice to save other aspects of the sampler state.  E.G., saving the state of the step-size adaptation process.  Maybe boost serialization would be a good option here after all.

[Bob Carpenter]

On 5/9/12 3:11 PM, Marcus Brubaker wrote:
> On Tue, May 8, 2012 at 8:52 PM, Andrew Gelman <gel...@stat.columbia.edu
> <mailto:gel...@stat.columbia.edu>> wrote:
>
> On May 8, 2012, at 8:25 PM, Marcus Brubaker wrote:
>
> > At the end of the day what constitutes a reasonable but diffuse
> random initialization is always going to be model (if not data)
> specific. Anything completely generic will be necessarily wrong most
> of the time :)
>
> Yes, but it's a good general principle to parameterize models so
> that parameters are of order 1.
>
>
> True, but sometimes that's hard to do unless you know a fair bit about
> the data in the first place.

Why not also normalize for rotation? Then we
can rescale the components independently :-)

> I recommend avoiding binaries! I'd prefer to have the object that
> is R-readable.

> As Bob noted, this would be only for restarting a chain that had
> finished or been interrupted. A binary blob isn't ideal but it is the
> easiest and fastest to implement. Now that I think about it, a more
> structured data format may be needed though, since it would be nice to
> save other aspects of the sampler state.

Good point. Though it'd be interesting to also allow
warmup/adaptation to be rerun for some amount of time.

In fact, that might even make sense in the ordinary
non-restarted course of things now that I think about it.
It will all depend on whether we can get to (and know we
get to) a reasonable state during warmup.

> E.G., saving the state of the
> step-size adaptation process. Maybe boost serialization would be a good
> option here after all.

If it's anything like Java serialization, I like it
already!

- Bob

[Marcus Brubaker]

I'm not sure how Java serialization works but boost serialization is pretty simple to use.  Here is an example from a different project I'm working on which is representing map data.

  struct Point {
    double lat,lon,alt;
    Point() {}
    Point(double lat, double lon, double alt) : lat(lat), lon(lon), alt(alt) {}
    template<typename A> void serialize(A &ar, const unsigned int){
      ar & BOOST_SERIALIZATION_NVP(lat)
         & BOOST_SERIALIZATION_NVP(lon)
         & BOOST_SERIALIZATION_NVP(alt);
    }
  };

This function handles both reading and writing data. and the interface goes something like:

std::ofstream ofs(filename);
boost::archive::binary_oarchive oa(ofs);
os << object;

You can also create text or xml archives.

The serialization is able to handle pointers and STL classes correctly and supports format versioning.  I'm not sure what it would do with Eigen classes, but I can't imagine that being insurmountable if it came up.

Cheers,
Marcus

[Bob Carpenter]

That sounds just like Java versioning other than
the funky use of the '&' operator and the macros, of course.
Presumably the template A is for the stream?

I'll read up on it later and maybe kick the tires.

I think we can get away with pretty simple serialization
of std::vector of primitive types.

- Bob

> <mailto:gel...@stat.columbia.__edu

[Marcus Brubaker]

Yeah, the template A is for either an input or output archive stream.

The & operator threw me for a loop initially too.  It is a stand-in for << or >> depending on whether A is an input or output stream.

Cheers,
Marcus