Stepsize Jitter Back

[Michael Betancourt]

I added step size jitter back into the HMC samplers, so we should have the exact same functionality as before.

High up on the to-do list is improved command line arguments, but that would require that we first (a) define the syntax for sub-arguments and then (b) update stan::io::cmd_line accordingly. Anyone have any opinions on (a)?

-Mike

[Jiqiang Guo]

There were two other arguments, delta and gamma.  Do we not need them now?

--

Jiqiang 

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

[Michael Betancourt]

Ah, yes.  Those are for the dual averaging.  We should either remove them or exposé all of the dual averaging parameters.

[Bob Carpenter]

We should leave those there until we want to
break backward compatibility.

Then we can expose more if anyone wants to.

Marcus had a good idea as to redefining syntax,
which is to have defaults with structure:

stan sample [general sampler args]
stan optimize [general optimizer args]

and then specific subcommands for users who want to
control the sampler parameters

stan hmc [hmc args] [general sampler args]
stan metropolis [metropolis args] [general sampler args]
...

stan bfgs [bfgs args] [general optimizer args]
...


BUT, we don't want to do this until Stan 2.0.0 if
it's at all possible to keep things backward
compatible for the next release.

IF we can sort out good options for all of this now,
then we could push up 2.0.0 to whenever the new
versions are completed.

Then 2.1.0 or 3.0.0 can be the RM-HMC release. Then people
will really be on our case about version creep!

- Bob

On 4/20/13 11:38 AM, Michael Betancourt wrote:
> Ah, yes. Those are for the dual averaging. We should either remove them or expos� all of the dual averaging parameters.

>
> On Apr 20, 2013, at 3:37 PM, Jiqiang Guo <guo...@gmail.com <mailto:guo...@gmail.com>> wrote:
>
>> There were two other arguments, delta and gamma. Do we not need them now?
>>
>> --
>> Jiqiang
>>
>>
>> On Sat, Apr 20, 2013 at 9:35 AM, Michael Betancourt <betan...@gmail.com <mailto:betan...@gmail.com>> wrote:
>>
>> I added step size jitter back into the HMC samplers, so we should have the exact same functionality as before.
>>
>> High up on the to-do list is improved command line arguments, but that would require that we first (a) define the
>> syntax for sub-arguments and then (b) update stan::io::cmd_line accordingly. Anyone have any opinions on (a)?
>>
>> -Mike
>>
>> --
>> 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>.

>> 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+u...@googlegroups.com>.

[Michael Betancourt]

> We should leave those there until we want to
> break backward compatibility

Okay, I'll implement the two options -- right now they're accepted but don't do anything.

> Marcus had a good idea as to redefining syntax,
> which is to have defaults with structure:
>
> stan sample [general sampler args]
> stan optimize [general optimizer args]
>
> and then specific subcommands for users who want to
> control the sampler parameters
>
> stan hmc [hmc args] [general sampler args]
> stan metropolis [metropolis args] [general sampler args]
> ...
>
> stan bfgs [bfgs args] [general optimizer args]

> …

Yup, totally agree. But how robust do we want it? Do we find some syntax (or assume an ordering) that associates subcommands with commands or just dump all commands and subcommands together? It seems like in either case we'll have to make significant changes to the stan::io::cmd_line, and I just wanted to settle on a format before beginning to look how to code it.

I have no opinion on versioning, but my first priority is diagnostics, then arguments, then RMHMC.

[Bob Carpenter]

Did the way the diagonal step sizes get estimated
change?

On 4/20/13 7:21 PM, Michael Betancourt wrote:
>> We should leave those there until we want to
>> break backward compatibility
>
> Okay, I'll implement the two options -- right now they're accepted but don't do anything.

???

I realize it's now a diagonal mass matrix instead
of varying step sizes (I'll have to change the manual
to reflect that).

But did the estimation strategy change so the dual
averaging params are no longer relevant? If so, it
doesn't make sense to leave them in place.

>> Marcus had a good idea as to redefining syntax,
>> which is to have defaults with structure:
>>
>> stan sample [general sampler args]
>> stan optimize [general optimizer args]
>>
>> and then specific subcommands for users who want to
>> control the sampler parameters
>>
>> stan hmc [hmc args] [general sampler args]
>> stan metropolis [metropolis args] [general sampler args]
>> ...
>>
>> stan bfgs [bfgs args] [general optimizer args]

>> �

>
> Yup, totally agree. But how robust do we want it? Do we find some syntax (or assume an ordering) that associates subcommands with commands or just dump all commands and subcommands together? It seems like in either case we'll have to make significant changes to the stan::io::cmd_line, and I just wanted to settle on a format before beginning to look how to code it.
>
> I have no opinion on versioning,

It's not so much versioning as how big a step we want
to take with each release. The version number follows
(with measurement error) from the scope of change.

> but my first priority is diagnostics, then arguments, then RMHMC.

And it won't really be "stan" from the command line,
but "my_model sample".

The executable is just "my_model", as now, but the first
thing it'll do is look at what the next argument is,
which would be required to be one of the subcommands.
That is, it should work like subcommands in git or svn.

It'd also be nice to have it check options, so that
unrecognized options are rejected.

Along those lines, the dual averaging params could be
accepted, but then spit out an error.

- Bob

[Michael Betancourt]

> I realize it's now a diagonal mass matrix instead
> of varying step sizes (I'll have to change the manual
> to reflect that).

Yes, when switching to a proper metric instead of the "varying step sizes" the step-size and metric elements became identifiable and so the adaptation had to be factored. Basically the variance (or covariance for dense) is estimated using the Welford algorithm while the step-size is estimated using dual-averaging.

> But did the estimation strategy change so the dual
> averaging params are no longer relevant? If so, it
> doesn't make sense to leave them in place.

No, they're used but take default values from internal constructors. It was just odd that we exposed two of the dual-averaging parameters but not the others, so I ignored the option to set them by hand when refactoring. I'd rather we make them all available, which the subcommands would make easy. For now I'll de-neuter the two existing options on master.

>
> And it won't really be "stan" from the command line,
> but "my_model sample".
>
> The executable is just "my_model", as now, but the first
> thing it'll do is look at what the next argument is,
> which would be required to be one of the subcommands.
> That is, it should work like subcommands in git or svn.
>
> It'd also be nice to have it check options, so that
> unrecognized options are rejected.

Yup, exactly what I was thinking. But how exactly do we implement it? Right now stan::io::cmd_line is agnostic to the commands, it just accumulates them and allows the user to query but any subcommand functionality would either require listing all of the possible command/subcommand relationships in stan::io::cmd_line or changing the syntax of the arguments. For example,

1) Denote subcommands with a prefix, such as s_

my_model sample hmc s_hmc_arg1 s_hmc_arg2 da s_da_arg1

This way, as stan::io::cmd_line reads the arguments it can associate commands (unprefixed) with subcommands (everything with the prefix immediately following the unprefixed command). But it's really ugly.

2) Denote subcommands with equals

my_model sample hmc:hmc_arg1,hmc_arg2 da:da_arg1

Again, easy to parse without knowing the command/subcommand relationships. Less ugly, but still pretty ugly.

3) Denote commands as flags

my_model sample -hmc hmc_arg1 hmc_arg2 -da da_arg1

4) Lookup table after parsing:

my_model sample hmc hmc_arg1 hmc_arg2 da da_arg1

gets parsed to a list

args=[hmc, hmc_arg1, hmc_arg2, da, da_arg1]

Then you just run over the list checking against a look-up table of commands and their subcommands.

5) Lookup table before parsing:

my_model sample hmc hmc_arg1 hmc_arg2 da da_arg1

While parsing you check against the big look-up table (which would have to be stored in stan::io::cmd_line now, unless we make a big giant class just to store the info). This way you can also restrict valid subcommands to those immediately following the command.

I'm just not sure what everyone's preference is. 3 seems to be the standard SVN approach although 5 wouldn't be that hard.

[Matt Hoffman]

I'm fine with whatever adaptation strategy works best. But we should
really have a way of determining that before making changes to core
functionality.

Also, I don't think the choice between multiple step sizes/diagonal
mass matrices in any way affects the identifiability question. The
position and energy trajectories are *in no way* changed by that
implementation choice. Any identifiability issues associated with one
implementation are every bit as present in the other; even if they're
obscured, they'll still affect the speed with which the adaptation
converges.

I do think it's important to factorize the questions of overall and
relative step size estimation, which means scaling the mass matrix so
that something (largest eigenvalue? smallest eigenvalue? trace?) is
constant.

Also, I think that even with the current strategy we're much too
sensitive to the early iterations when fitting the per-variable step
sizes. Estimates of the covariance matrix from the initial warmup
samples are garbage—they only tell you how far away you initialized
the sampler from a mode. I'd like to have the mass matrix adaptation
turned off until we have some reason to believe that things have kind
of stabilized.

Matt

[Michael Betancourt]

> Also, I don't think the choice between multiple step sizes/diagonal
> mass matrices in any way affects the identifiability question. The
> position and energy trajectories are *in no way* changed by that
> implementation choice. Any identifiability issues associated with one
> implementation are every bit as present in the other; even if they're
> obscured, they'll still affect the speed with which the adaptation
> converges.
>
> I do think it's important to factorize the questions of overall and
> relative step size estimation, which means scaling the mass matrix so
> that something (largest eigenvalue? smallest eigenvalue? trace?) is
> constant.

This is not exactly true. The problem is that, in general, you can't scale the metric arbitrarily; only in the specific case of the explicit leapfrog integrator do you have that freedom. In that case you can perform a symplectomorphism

P = (L^{T})^{-1} p
Q = L^{-1} q

and then write the integrator in terms of P and q so that you always have (epsilon * (L^{T})^{-1}) grouped together, and you're free to exchange scales arbitrarily. Because the refactoring was meant to separate the sampling algorithm, the Hamiltonian, the integrator, and the potential such a strategy was not feasible. In particular, it wouldn't work for RMHMC and implementing it would require an architecture necessitating rewrites for each new algorithm.

That said, the only real difference in the adaptation is that an update in the metric now necessitates a corresponding update of the step size or else the acceptance rate will be off. Given the performance of your step size initialization routine, however, that's hardly a problem.

> Also, I think that even with the current strategy we're much too
> sensitive to the early iterations when fitting the per-variable step
> sizes. Estimates of the covariance matrix from the initial warmup
> samples are garbage—they only tell you how far away you initialized
> the sampler from a mode. I'd like to have the mass matrix adaptation
> turned off until we have some reason to believe that things have kind
> of stabilized.

Agreed. Implementing a better warm-up is on my very long, and rapidly growing, to-do list. Hopefully once I'm settled I'll be able to burn through a bunch of it quickly and get to the warm-up sooner rather than later.

[Matt Hoffman]

I see. Sounds like it should be possible to implement it in a way that
hides the effect of overall mass matrix scaling from the step size
adaptation algorithm, so that should be fine.

[Bob Carpenter]

On 4/20/13 8:20 PM, Michael Betancourt wrote:
...

> Basically the variance (or covariance for dense) is estimated using
> the Welford algorithm while the step-size is estimated using dual-averaging.

...

Would you mind breaking out the accumulator into its
own class for testing and re-use purposes?

Something like:

class welford_accumulator {
private:
... data to store ...
public:
welford_accumulator(double prior_count);
void add_sample(const vector& y);
long num_samples() const;
void sample_mean(vector& mu) const;
void sample_variance(vector& sigma_sq) const;
void sample_covariance(matrix& Sigma) const;
};

Can the regularization you're doing,

(N*Sigma + 5*I) / (N+5),

be coded up more efficiently by just initialzing the counters?
That is, treating it like prior data?


That's why I called 5 the prior count. I also don't know
technically if that should also be included in the N * (N-1)
adjustment for bias term -- looks like your code didn't
include it.

I think Yuanjun was using the estimator

(N*Sigma + e1*diag(Sigma) + e2*I) / (N + e1 + e2)

Any reason for the switch?

This accumulator class could be templated on the real type,
but I think we'd probably only ever use the double version.
I don't see plugging our autodif vars into it.

- Bob

[Michael Betancourt]

> Would you mind breaking out the accumulator into its
> own class for testing and re-use purposes?
>
> Something like:
>
> class welford_accumulator {
> private:
> ... data to store ...
> public:
> welford_accumulator(double prior_count);
> void add_sample(const vector& y);
> long num_samples() const;
> void sample_mean(vector& mu) const;
> void sample_variance(vector& sigma_sq) const;
> void sample_covariance(matrix& Sigma) const;
> };

It is already. var_adapter does variance accumulation/estimation for a vector of parameters while covar_adapter does covariance accumulation/estimation, and that's it. They have no other purpose. They'll fall right into Metropolis adaptation and the like. How else would you factor things?

> Can the regularization you're doing,
>
> (N*Sigma + 5*I) / (N+5),
>
> be coded up more efficiently by just initialzing the counters?
> That is, treating it like prior data?

> That's why I called 5 the prior count. I also don't know
> technically if that should also be included in the N * (N-1)
> adjustment for bias term -- looks like your code didn't
> include it.

I'm not sure if the Welford estimation commutes with regularization, in any case any speed up would be pretty miniscule.

> I think Yuanjun was using the estimator
>
> (N*Sigma + e1*diag(Sigma) + e2*I) / (N + e1 + e2)
>
> Any reason for the switch?

Just to unify the variance/covariance implementations while I was coding everything up. It's really easy to change.

[Michael Betancourt]

> And it won't really be "stan" from the command line,
> but "my_model sample".
>
> The executable is just "my_model", as now, but the first
> thing it'll do is look at what the next argument is,
> which would be required to be one of the subcommands.
> That is, it should work like subcommands in git or svn.
>
> It'd also be nice to have it check options, so that
> unrecognized options are rejected.
>
> Along those lines, the dual averaging params could be
> accepted, but then spit out an error.

I mocked up what I thought could work. You can try it out be checking out the arg_refactor branch and building any model (nuts_command exits before any actual sampling code gets called). Alternatively, the new code is all in src/stan/gm/arguments and is called a la

std::vector<argument*> valid_arguments;
valid_arguments.push_back(new arg_static_hmc());

argument_reader reader(valid_arguments);

reader.parse_arguments(argc, argv);
reader.print_arguments(&std::cout);

Example output:

[~/Documents/Research/Code/Inference/stan-dev/stan (arg_refactor)]-betan: ./example -static_hmc stepsize=1.0 old_nuts int_time=-1.0
Found the sub-argument old_nuts without any value, it will be ignored.
WARNING: -1 is not a valid value for int_time, which must be positive
falling back to default value

int_time 0 (Default)
stepsize 1

Any thoughts?

[Bob Carpenter]

On 4/21/13 3:28 PM, Michael Betancourt wrote:
>> Would you mind breaking out the accumulator into its
>> own class for testing and re-use purposes?
>>
>> Something like:
>>
>> class welford_accumulator {
>> private:
>> ... data to store ...
>> public:
>> welford_accumulator(double prior_count);
>> void add_sample(const vector& y);
>> long num_samples() const;
>> void sample_mean(vector& mu) const;
>> void sample_variance(vector& sigma_sq) const;
>> void sample_covariance(matrix& Sigma) const;
>> };
>
> It is already. var_adapter does variance accumulation/estimation for a vector of parameters while covar_adapter does covariance accumulation/estimation, and that's it. They have no other purpose. They'll fall right into Metropolis adaptation and the like. How else would you factor things?

Close, but the online covariance estimator isn't fully
abstracted.

1. It extends the non-trivial base class mcmc::stepsize_adapter,
which isn't relevant for a covariance estimator.

2. It also calls stepsize_adapter::init() in its init()
function and the dependency isn't relevant.

3. There's a branch on whether to do adaptation or not,
which also isn't relevant.

4. It's in the wrong namespace (stan::mcmc). for a general

What I was suggesting is pulling out just the general
functionality as outlined in the public interface above.

You would keep your covar_adapter class with the same
functionality as it has now in mcmc. But it would hold
an instance of an online covariance estimator rather than
many of the individual variables it has now.

Also, why use a separate init() function? That's what the
constructor's for and it automatically calls the parent
class init() so you don't have to code it yourself. It
also lets you call an alternative constructor if you want.

>> Can the regularization you're doing,
>>
>> (N*Sigma + 5*I) / (N+5),
>>
>> be coded up more efficiently by just initialzing the counters?
>> That is, treating it like prior data?
>
>> That's why I called 5 the prior count. I also don't know
>> technically if that should also be included in the N * (N-1)
>> adjustment for bias term -- looks like your code didn't
>> include it.
>
> I'm not sure if the Welford estimation commutes with regularization, in any case any speed up would be pretty miniscule.

I wasn't thinking in terms of speed, but rather
in terms of code simplicity. Though the pseudo-data
approach might arguably be more confusing to someone
reading to the code -- how you have it makes it very
clear what's going on for estimation.

>> I think Yuanjun was using the estimator
>>
>> (N*Sigma + e1*diag(Sigma) + e2*I) / (N + e1 + e2)
>>
>> Any reason for the switch?
>
> Just to unify the variance/covariance implementations while I was coding everything up. It's really easy to change.

We should talk about this with Yuanjun and Andrew
at the meeting. I suspect we don't have enough comparisons
to come to any sound conclusions yet --- benchmarking with
the noise introduced by MCMC plus random inits is a challenge.

- Bob

[Bob Carpenter]

I'd prefer to halt if the command is misspecified
rather than guessing the user's intent and continuing.

Opinions on this vary widely, I think. Anyone else
have an opinion?

We should think hard about the arg names. The current
ones aren't great. What is int_time?

I don't like the "negative number means infinity" kind
of things we're doing now (and I did them!). I'm not
sure what a good alternative idiom is. Maybe a different
command? That is, I don't like this choice:

--epsilon=0 // adapt eps
--epsilon=.01 // use eps=.01

I would prefer:

--adapt_epsilon // adapt
--epsilon=.01 // use eps = .01
--adapt_epsilon --epsilon=.01 // ILLEGAL!

Thoughts?

I really don't have a lot of experience in designing command-line
or R-style interfaces.

For defaults, can we have max defaults be the the max
legal value rather than some value meaning no max?

I prefer "./example static_hmc..." to "./example -static_hmc".

I like the double dashes used by git for command property args,
but am not wedded to it. I did it originally partly because
it seemed easier to code.

- Bob

[Michael Betancourt]

> I'd prefer to halt if the command is misspecified
> rather than guessing the user's intent and continuing.
>
> Opinions on this vary widely, I think. Anyone else
> have an opinion?

I have no problem either way, provided that in either case the user gets a clear error message so they know what to fix.

> We should think hard about the arg names. The current
> ones aren't great. What is int_time?

Total integration time for static HMC (I'm going to be pretty insistent on step size and integration time begin the free parameters for static HMC instead of step size and number of steps). But in any case the names are easy fixes.

> I don't like the "negative number means infinity" kind
> of things we're doing now (and I did them!). I'm not
> sure what a good alternative idiom is. Maybe a different
> command? That is, I don't like this choice:
>
> --epsilon=0 // adapt eps
> --epsilon=.01 // use eps=.01
>
> I would prefer:
>
> --adapt_epsilon // adapt
> --epsilon=.01 // use eps = .01
> --adapt_epsilon --epsilon=.01 // ILLEGAL!
>
> Thoughts?

Agreed! I think we can fully factor the various options so that we don't have to overload arguments like this. For example,

./example sample --iter=2000 --warmup=500 nuts --max_depth=10 --stepsize=0.1 --metric=diag_e adapt --delta=0.65 output --sample_file=samples.csv --diagnostic_file=diagnostic.csv

Then in the parsing we check to determine the final algorithm. NUTS (check). Diagonal Euclidean metric (check). Adaptation (check) -> adapt_diag_e_nuts. This is what I had in mind when building the mock-up.

> For defaults, can we have max defaults be the the max
> legal value rather than some value meaning no max?

I think we should be careful about that. For example, in the case of max_depth for NUTS we don't want misbehaving models to build an infinitely large tree.

> I prefer "./example static_hmc..." to "./example -static_hmc".
>
> I like the double dashes used by git for command property args,
> but am not wedded to it. I did it originally partly because
> it seemed easier to code.

Trivial change (see above for example). Again, I have no strong preference other than not building code that people don't like because dashes are in the wrong place.

[Michael Betancourt]

> Close, but the online covariance estimator isn't fully
> abstracted.
>
> 1. It extends the non-trivial base class mcmc::stepsize_adapter,
> which isn't relevant for a covariance estimator.
>
> 2. It also calls stepsize_adapter::init() in its init()
> function and the dependency isn't relevant.
>
> 3. There's a branch on whether to do adaptation or not,
> which also isn't relevant.
>
> 4. It's in the wrong namespace (stan::mcmc). for a general
>
> What I was suggesting is pulling out just the general
> functionality as outlined in the public interface above.
>
> You would keep your covar_adapter class with the same
> functionality as it has now in mcmc. But it would hold
> an instance of an online covariance estimator rather than
> many of the individual variables it has now.

Ah, I see.

(1, 2) Easily fixed by just having the adaptive algorithms derive from both stepsize_adapter and (co)var_adapter. Probably something I should do anyways.
(4) What namespace would be considered appropriate? math?

I can do it, just so long as things don't go crazy and someone makes me abstract the regularization away, too. At least not until we use it more than once.

> Also, why use a separate init() function? That's what the
> constructor's for and it automatically calls the parent
> class init() so you don't have to code it yourself. It
> also lets you call an alternative constructor if you want.

It's there in case you want to restart the adaptation, with multi-stage warmup algorithms in mind. I only use it in the constructor since it's there anyways.

> I wasn't thinking in terms of speed, but rather
> in terms of code simplicity. Though the pseudo-data
> approach might arguably be more confusing to someone
> reading to the code -- how you have it makes it very
> clear what's going on for estimation.

Right; that was my focus in cleaning up the algorithms.

> We should talk about this with Yuanjun and Andrew
> at the meeting. I suspect we don't have enough comparisons
> to come to any sound conclusions yet --- benchmarking with
> the noise introduced by MCMC plus random inits is a challenge.

Sure. I hope that the refactoring makes comparisons really easy.

[Bob Carpenter]

I think stan::prob makes more sense than stan::math.

I'd rather not introduce multiple inheritance of
implementations if we don't have to.

I'd much rather see you just include an instance
(not a reference to an instance) as a member variable.

As a design strategy, containing is much more flexible
than inheriting because you can choose what to expose.
It does mean you need to define the
methods again where you just need to pass-through
results from the accumulator. But you probably want
to name them something different anyway.

If you do multistage warmups, you probably do NOT
want to call init() twice. Or you should at least
name it something other than init(). Maybe an
init() to do allocation, etc. and a reset()?
Then init() can call reset()?

With lots of structures in programs, it's much easier
to reason about them if they get constructed once.
It's not that expensive to construct a new object
and a lot cleaner than tweaking an existing one.
The tweak needs to be sure to clear all the state,
and if done incrementally, to insure consistency.

- Bob

[Bob Carpenter]

On 4/21/13 6:18 PM, Michael Betancourt wrote:
>> I'd prefer to halt if the command is misspecified
>> rather than guessing the user's intent and continuing.
>>
>> Opinions on this vary widely, I think. Anyone else
>> have an opinion?
>
> I have no problem either way, provided that in either case the user gets a clear error message so they know what to fix.

Absolutely the top priority!

>> We should think hard about the arg names. The current
>> ones aren't great. What is int_time?
>
> Total integration time for static HMC (I'm going to be pretty insistent on step size
> and integration time begin the free parameters for static HMC instead of
> step size and number of steps). But in any case the names are easy fixes.

:-) Makes sense to me. We'll have to explain this pretty
clearly. Static HMC is the alternative to NUTS HMC?

Are we going to name something to do adaptation during
warmup?

Are you changing all of these things now or are we
talking about a subsequent release?

I can't even figure out whether sooner or later is better
for RStan, PyStan, and our poor users who are getting
a lot thrown at them.

>> I don't like the "negative number means infinity" kind
>> of things we're doing now (and I did them!). I'm not
>> sure what a good alternative idiom is. Maybe a different
>> command? That is, I don't like this choice:
>>
>> --epsilon=0 // adapt eps
>> --epsilon=.01 // use eps=.01
>>
>> I would prefer:
>>
>> --adapt_epsilon // adapt
>> --epsilon=.01 // use eps = .01
>> --adapt_epsilon --epsilon=.01 // ILLEGAL!
>>
>> Thoughts?
>
> Agreed! I think we can fully factor the various options so that we don't have to overload arguments like this. For example,
>
> ./example sample --iter=2000 --warmup=500 nuts --max_depth=10 --stepsize=0.1 --metric=diag_e adapt --delta=0.65 output --sample_file=samples.csv --diagnostic_file=diagnostic.csv
>
> Then in the parsing we check to determine the final algorithm. NUTS (check). Diagonal Euclidean metric (check). Adaptation (check) -> adapt_diag_e_nuts. This is what I had in mind when building the mock-up.
>
>> For defaults, can we have max defaults be the the max
>> legal value rather than some value meaning no max?
>
> I think we should be careful about that. For example, in the case of max_depth for NUTS we don't want misbehaving models to build an infinitely large tree.

Oh, I just meant defaults that were equal to some
large value. For NUTS, we should be careful. And
probably also print a warning to the output stream
when it hits the max depth as a diagnostic.

>> I prefer "./example static_hmc..." to "./example -static_hmc".
>>
>> I like the double dashes used by git for command property args,
>> but am not wedded to it. I did it originally partly because
>> it seemed easier to code.
>
> Trivial change (see above for example). Again, I have no strong preference other than not building code that people don't like because dashes are in the wrong place.

I don't have even a weak preference. I used the
double dashes because it was easy to code.

- Bob

[Michael Betancourt]

> :-) Makes sense to me. We'll have to explain this pretty
> clearly. Static HMC is the alternative to NUTS HMC?

Yessir.

> Are we going to name something to do adaptation during
> warmup?

See the example -- my though is to have adaptation switched with its own flag (adapt=on --delta=0.65). Figuring out the exact algorithm to use would be a post-argument processing step.

> Are you changing all of these things now or are we
> talking about a subsequent release?

The sampler functionality is there, just waiting on better arguments. So not until we've redone the arguments.

> I can't even figure out whether sooner or later is better
> for RStan, PyStan, and our poor users who are getting
> a lot thrown at them.

Understandable, but I'm focusing on a schema that enables automated documentation which, along with the verbose error handling, should help users.

> Oh, I just meant defaults that were equal to some
> large value. For NUTS, we should be careful. And
> probably also print a warning to the output stream
> when it hits the max depth as a diagnostic.

Sorry -- misread your previous comment as "max floating point value" not "well-chosen maximum". :-)

> I don't have even a weak preference. I used the
> double dashes because it was easy to code.

I was going over the code last night and convinced myself that, to really make this robust, we'd want to implement recursive argument handling (hmc is a subargument of sample, but hmc has its own options like algorithm=nuts, which has its own options). I've sketched out some code that should enable this functionality but such deep hierarchies might not be clear on the command line. So does anyone have any preference between

(flat) ./example method=sample adapt=on delta=0.65 algorithm=hmc metric=diag_e variant=nuts max_depth=10

verses

(hierarchical flags) ./example method=sample -adapt=on --delta=0.65 -algorithm=hmc --metric=diag_e --variant=nuts ---max_depth=10

?

[Michael Betancourt]

On Apr 22, 2013, at 2:18 AM, Bob Carpenter <ca...@alias-i.com> wrote:

> I think stan::prob makes more sense than stan::math.
>
> I'd rather not introduce multiple inheritance of
> implementations if we don't have to.
>
> I'd much rather see you just include an instance
> (not a reference to an instance) as a member variable.
>
> As a design strategy, containing is much more flexible
> than inheriting because you can choose what to expose.
> It does mean you need to define the
> methods again where you just need to pass-through
> results from the accumulator. But you probably want
> to name them something different anyway.

Oh, I was referring to multiple inheritance in the adaptive samplers. For example

class adapt_diag_e_nuts: public diag_e_nuts, public stepsize_adapter, public var_adapter

so that var_adapter doesn't have to inherent from stepsize_adapter unnecessarily. var_adapter would hold the instance of welford_var_estimator

> If you do multistage warmups, you probably do NOT
> want to call init() twice. Or you should at least
> name it something other than init(). Maybe an
> init() to do allocation, etc. and a reset()?
> Then init() can call reset()?
>
> With lots of structures in programs, it's much easier
> to reason about them if they get constructed once.
> It's not that expensive to construct a new object
> and a lot cleaner than tweaking an existing one.
> The tweak needs to be sure to clear all the state,
> and if done incrementally, to insure consistency.

Sure. The naming referred to an algorithmic initialization, not an instance initialization. But I can see the confusion and will make the name more explicit.

[Bob Carpenter]

On 4/22/13 4:28 AM, Michael Betancourt wrote:
>
> On Apr 22, 2013, at 2:18 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>

...

>> As a design strategy, containing is much more flexible
>> than inheriting because you can choose what to expose.
>> It does mean you need to define the
>> methods again where you just need to pass-through
>> results from the accumulator. But you probably want
>> to name them something different anyway.
>
> Oh, I was referring to multiple inheritance in the adaptive samplers. For example
>
> class adapt_diag_e_nuts: public diag_e_nuts, public stepsize_adapter, public var_adapter
>
> so that var_adapter doesn't have to inherent from stepsize_adapter unnecessarily. var_adapter would hold the instance of welford_var_estimator

That makes much more sense to me because these are all
is-a relations.

I'd still rather not do multiple inheritance of implementations
unless we absolutely have to. For one, I'm not even sure how
it works! But I can learn...

- Bob

[Michael Betancourt]

Adaptation refactored (Welford estimator classes, adaptation classes, adapter classes that wrap the adaptation algorithms), tested, and pushed.