Re: Long-ish email with questions about Stan

[Bob Carpenter]

[moved to stan-dev]

On Jul 22, 2014, at 9:48 PM, Alp Kucukelbir <a...@cs.columbia.edu> wrote:

> I'm making progress! Faster than I had thought, actually.

Great news.

> I have three (direct) questions for you guys and one favor to ask:
>
> [Question 1:] When I do
>
> stan::model::gradient(model_, cont_params_, log_prob, init_grad, &std::cout);
>
> versus
>
> log_prob = model_.template log_prob<false, true>(cont_params_, &std::cout);
>
> I get slightly different (yet close) numerical values for log_prob.
>
> Is there a numerical difference between the two methods? Or am I missing the Jacobian somewhere?

The three template parameters for log_prob are:

propto (boolean) whether to drop constant terms
jacobian: (boolean) whether to include the Jacobian
T: (scalar) type of scalars (stan::agrad::var, double, etc.)

The last one will be inferred automatically based on the argument

The gradient method sets both booleans to true.

We need to get the model concept better documented. The current
doc is here:

https://github.com/stan-dev/stan/wiki/Model-Concept

Please update as you find errors. For instance, I just noticed the
scalar type parameter T is in the wrong order on log_prob (so
I fixed it); there's probably a lot more that's not up to date
or correct in it.

> [Question 2:] I know that I can call stan::model::gradient to get the gradient of the log prob. The Jacobian is taken care of as well.
>
> But how should I go about calculating the gradient of the log prob composed with some affine function g(params)? I have the function implemented in C++, but I suspect I'll have to do some agrad magic...
>
> Daniel had sent me this elucidating example (see attached). Should I follow that? Or is there some built in functionality for what I'm trying to accomplish already built into Stan?

Assuming I'm understanding correctly and you want to calculate

d/dparams log_prob(g(params)),

there are a bunch of ways you can do it.

I'd suggest using the functional in stan/agrad/autodiff.hpp.

Just define a new functor class

template <typename M>
struct affine_log_prob {
const M& model_;

affine_log_prob(const M& model) : model_(model) { }

template <typename T>
T operator()(Matrix<T,Dynamic,1> x) {
return model_.template log_prob<true,true>(g(params));
}
};

Then just pass an instance of affine_log_prob into the gradient
function. It's important to keep the reference to model_ so you
don't inadvertently make a copy --- model instances have all the data
instantiated.

You'll need to construct the instance with both the model (and specify
the appropriate type to the struct's template parameter) and
also whatever parameters are involved in computing g (e.g., a matrix
and vector if it's an affine transform with translation).

> [Question 3:] While I'm developing on /feature/bbvb/, I would like to "keep tabs" on what's happening on /develop. Currently I do the following
>
> (I'm working on /feature/bbvb.)
> git pull
> git merge develop
>
> But I read that I might want to be doing "rebasing" instead. What do you guys recommend?

Rebasing erases the history of your commits, which is something
different. We want to make sure all the commits at least compile so
that we can run git bisect to find problems when necessary. If they
don't, then rebasing can be used --- I don't even know how to do it.

You want to do:

> git pull origin develop

in order to merge the latest develop into your feature.

> [Favor 1:] Could one of you take a look at /feature/bbvb/ and tell me how I'm ruining all that is beautiful in Stan? I would appreciate any tips/pointers/suggestions/etc on the code.

What should I look at in that branch (feature/bbvb --- you
don't need leading/trailing slashes)?

The main thing that code needs is to not be in too quirky a style
and to pass reasonable unit tests.

- Bob

[Alp Kucukelbir]

> The three template parameters for log_prob are:
>
> propto (boolean) whether to drop constant terms
>
> jacobian: (boolean) whether to include the Jacobian
>
> T: (scalar) type of scalars (stan::agrad::var, double, etc.)

Yes, but even when I match them, I get a small numerical difference. E.g.

log_prob =
-11251.1

log_prob =
-11710.8

> We need to get the model concept better documented. The current
> doc is here:
> https://github.com/stan-dev/stan/wiki/Model-Concept

Thanks for this pointer.

> Please update as you find errors.

I'll try my best. (I will ask you guys on Thursday about how to go about fixing small things like this: branches, hotfixes, etc.)

> Assuming I'm understanding correctly and you want to calculate
> d/dparams log_prob(g(params)),
> there are a bunch of ways you can do it.
>
> I'd suggest using the functional in stan/agrad/autodiff.hpp.

OK, thanks for this! Michael also suggested just hardcoding the derivative of the affine transform, since it's so simple. But I'm not sure if we're only going to stick with affine transforms... so for anything more complicated, I'll take the approach you suggest.

> Rebasing erases the history of your commits, which is something
> different. We want to make sure all the commits at least compile so
> that we can run git bisect to find problems when necessary. If they
> don't, then rebasing can be used --- I don't even know how to do it.

OK. I don't fully understand, but I'll stay away from rebasing for now.

> What should I look at in that branch (feature/bbvb --- you
> don't need leading/trailing slashes)?

Mainly the stuff I'm adding to /src/stan/vb and modifications to /src/stan/common/command.hpp

> The main thing that code needs is to not be in too quirky a style
> and to pass reasonable unit tests.

OK great.

Thanks guys!
Alp

[Bob Carpenter]

On Jul 23, 2014, at 4:05 AM, Michael Betancourt <betan...@gmail.com> wrote:

>> [Question 1:] When I do
>>
>> stan::model::gradient(model_, cont_params_, log_prob, init_grad, &std::cout);
>>
>> versus
>>
>> log_prob = model_.template log_prob<false, true>(cont_params_, &std::cout);
>>
>> I get slightly different (yet close) numerical values for log_prob.
>>
>> Is there a numerical difference between the two methods? Or am I missing the Jacobian somewhere?
>

> To clarify Bob’s response, gradient calls
>
> log_prob = model_.template log_prob<true, true>(cont_params_, &std::cout);
>
> which adds the normalization constants of the distributions.
>

Michael's right about what the code's doing, but this isn't how
it should be called because we don't need to compute the normalizing
constant.

The model_functional class (which should be called "model_functor")
is indeed specifying log_prob<true,true>, but it should be calling
log_prob<false,true>. At least for HMC, NUTS and optimization.
So one option if you need it to be <true,true> elsewhere is to
add two more template parameters to model_functor (assuming it
gets renamed properly).

I'm creating an issue that I think we should deal with for 2.4.0++.

https://github.com/stan-dev/stan/issues/797

>> [Question 2:] I know that I can call stan::model::gradient to get the gradient of the log prob. The Jacobian is taken care of as well.
>>
>> But how should I go about calculating the gradient of the log prob composed with some affine function g(params)? I have the function implemented in C++, but I suspect I'll have to do some agrad magic...
>>
>> Daniel had sent me this elucidating example (see attached). Should I follow that? Or is there some built in functionality for what I'm trying to accomplish already built into Stan?
>

> If it’s affine then you should be able to modify the gradient directly,
>
> g = a * f + b
> dg/dx = a * df/dx + b,
>
> no? In any case, if g is nonlinear then you’ll have to implement a functor.

Nice to have someone who actually understands calculus!
Michael's suggestion should be more efficient than calculating
via autodiff.

>> [Question 3:] While I'm developing on /feature/bbvb/, I would like to "keep tabs" on what's happening on /develop. Currently I do the following
>>
>> (I'm working on /feature/bbvb.)
>> git pull
>> git merge develop
>>
>> But I read that I might want to be doing "rebasing" instead. What do you guys recommend?
>

> git pull -ff origin develop to keep up with dev. You would rebase your branch before merge back into develop.

You should only rebase if the intermediate commits do not leave
Stan in a coherent state (and you probably shouldn't be doing that).

If anyone can explain how to do rebasing in a safe way, could they
put it in the developer process wiki on stan-dev/stan? I don't even
know how to do it and it looks scary.

- Bob

[Bob Carpenter]

On Jul 23, 2014, at 8:50 AM, Alp Kucukelbir <akucu...@gmail.com> wrote:

>> The three template parameters for log_prob are:
>>
>> propto (boolean) whether to drop constant terms
>>
>> jacobian: (boolean) whether to include the Jacobian
>>
>> T: (scalar) type of scalars (stan::agrad::var, double, etc.)
>
> Yes, but even when I match them, I get a small numerical difference. E.g.
>
> log_prob =
> -11251.1
>
> log_prob =
> -11710.8

Could you send the code that's calling things this way. It's
impossible to diagnose with just numbers. Or feel free to point
me to a branch.

>> We need to get the model concept better documented. The current
>> doc is here:
>> https://github.com/stan-dev/stan/wiki/Model-Concept
>
> Thanks for this pointer.
>
>> Please update as you find errors.
>
> I'll try my best. (I will ask you guys on Thursday about how to go about fixing small things like this: branches, hotfixes, etc.)

You can just edit the wiki directly.

Everything else should be outlined in the developer process.
We've never used hotfixes because we keep develop in a releasable
state.

>> Assuming I'm understanding correctly and you want to calculate
>> d/dparams log_prob(g(params)),
>> there are a bunch of ways you can do it.
>>
>> I'd suggest using the functional in stan/agrad/autodiff.hpp.
>
> OK, thanks for this! Michael also suggested just hardcoding the derivative of the affine transform, since it's so simple. But I'm not sure if we're only going to stick with affine transforms... so for anything more complicated, I'll take the approach you suggest.

Sounds good.

>> Rebasing erases the history of your commits, which is something
>> different. We want to make sure all the commits at least compile so
>> that we can run git bisect to find problems when necessary. If they
>> don't, then rebasing can be used --- I don't even know how to do it.
>
> OK. I don't fully understand, but I'll stay away from rebasing for now.

Then make sure every commit (not every push) leaves Stan in
a buildable state.

>
>> What should I look at in that branch (feature/bbvb --- you
>> don't need leading/trailing slashes)?
>
> Mainly the stuff I'm adding to /src/stan/vb and modifications to /src/stan/common/command.hpp

Michael's probably the best one to look at this. I haven't played
with the command structure much and barely understand VB.

- Bob

[Alp Kucukelbir]

>> Yes, but even when I match them, I get a small numerical difference. E.g.
>>
>> log_prob =
>> -11251.1
>>
>> log_prob =
>> -11710.8
>
> Could you send the code that's calling things this way. It's
> impossible to diagnose with just numbers. Or feel free to point
> me to a branch.

Sure. Here's what I did:

- check out /feature/bbvb/
- use cmdstan to compile some model (e.g. bayes lin regression with ARD)
- call ./executable vb output file=bla data file=bla.data.R

This will execute the function test() in src/stan/vb/bbvb.hpp

Lines 72 to 105 have the calculation for both "versions" of the log_prob.

[Michael Betancourt]

Alp — checkout https://github.com/stan-dev/stan/pull/799. If you want to merge the changes then go ahead and press merge and delete the bbvb_tweak branch. Merging a branch into develop is a whole ‘nother story, but we’re far away from that.

> --
> 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/d/optout.

[Alp Kucukelbir]

Hi Michael,

Thanks so much. I'll take a look at these tonight and try to merge the changes.

I'll follow up with any questions/comments tomorrow during the meeting.

(Thanks for the gentle introduction to pull requests as well!)

Cheers,
Alp


On Jul 23, 2014, at 4:41 PM, Michael Betancourt <betan...@gmail.com> wrote:

> Alp -- checkout https://github.com/stan-dev/stan/pull/799. If you want to merge the changes then go ahead and press merge and delete the bbvb_tweak branch. Merging a branch into develop is a whole 'nother story, but we're far away from that.

> You received this message because you are subscribed to a topic in the Google Groups "stan development mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-dev/KaKki-sX6Uk/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-dev+u...@googlegroups.com.

[Michael Betancourt]

Alp,

With what model file are you using to complete the gradients?

[Alp Kucukelbir]

Bayesian Linear Regression with ARD.

I ran this model by Daniel a while back. He fixed a bunch of problems with my initial version.

Seems to work with sampling and optimization. I thought it would be a good testbed for VB as the model has no local latent variables, but exhibits a relatively complex global latent variable structure.

bayes_lin_reg_ard.stan: stan model
bayes_lin_reg_ard.src.R: R code to generate data and "true parameters" files

bayes_lin_reg_ard.data.R: file with data to be loaded into Stan
bayes_lin_reg_ard.truth.R: file with "true" configuration of parameters (for comparing Stan's output)

[Michael Betancourt]

Alp,

Firstly, your calls to model::log_prob are inconsistent. You call model::log_prob<false, true> and model::log_prob<true, false>, but not
model::log_prob<true, true> which what gradient calls. Now you will still see a discrepancy between the log_prob you get from
gradient and the one you get from model::log_prob<true, true> because the latter treats everything as a constant, so only the
Jacobian survives. In the former only terms that don’t depend on the inputs are dropped.

Once you take this into account everything is consistent.

> <bayes_lin_reg_ard.data.R><bayes_lin_reg_ard.src.R><bayes_lin_reg_ard.stan><bayes_lin_reg_ard.truth.R>--

[Andrew Gelman]

“ARD”?? Let’s try to use our multilevel modeling jargon…!

> <bayes_lin_reg_ard.data.R><bayes_lin_reg_ard.src.R><bayes_lin_reg_ard.stan><bayes_lin_reg_ard.truth.R>--

[Alp Kucukelbir]

Firstly: my apologies, Andrew, for this (and the inevitable) future jargon offenses. I am but a mere image processing guy masquerading as a statistician. :)

Secondly: re: the log prob business, is the consensus that I use these two calls (when necessary)?

stan::model::gradient(model_, cont_params_, lp, init_grad, &std::cout);
lp = model_.template log_prob<true, true>(cont_params_, &std::cout);

Finally: re being a good citizen, not setting the world ablaze, etc... what is it exactly that I should be doing as I implement variational bayes?

[Bob Carpenter]

On Jul 25, 2014, at 3:16 PM, Alp Kucukelbir <akucu...@gmail.com> wrote:

> Firstly: my apologies, Andrew, for this (and the inevitable) future jargon offenses. I am but a mere image processing guy masquerading as a statistician. :)

I'm just a simple country computer scientist myself.

In this case, I don't know what Andrew's objecting to or
suggesting as an alternative --- Andrew, it'd help if you
were more concrete in what you want to see rather than just
name calling someone else's language as "jargon", which as
we know, is totally in the eye of the beholder (belistener?).

> Secondly: re: the log prob business, is the consensus that I use these two calls (when necessary)?

>
> stan::model::gradient(model_, cont_params_, lp, init_grad, &std::cout);
> lp = model_.template log_prob<true, true>(cont_params_, &std::cout);

Depends what you want to do. These will drop constant terms and apply
the Jacobian.

> Finally: re being a good citizen, not setting the world ablaze, etc... what is it exactly that I should be doing as I implement variational bayes?

I'm not sure what you're asking. Did you see the Developer Process
wiki page on stan-dev/stan? Follow that. We're happy to clarify if
it's unclear.

- Bob

[Andrew Gelman]

On Jul 25, 2014, at 9:23 PM, Bob Carpenter <ca...@alias-i.com> wrote:

>
> On Jul 25, 2014, at 3:16 PM, Alp Kucukelbir <akucu...@gmail.com> wrote:
>
>> Firstly: my apologies, Andrew, for this (and the inevitable) future jargon offenses. I am but a mere image processing guy masquerading as a statistician. :)
>
> I'm just a simple country computer scientist myself.
>
> In this case, I don't know what Andrew's objecting to or
> suggesting as an alternative

I think that “ARD” is what we call hierarchical or multilevel modeling in statistics. I have no objection to referring to the model as ARD also, I just don’t want ARD to be the main title. I think ARD stands for “automatic relevance determination” but I don’t really like that expression, as I don’t really see it as any more automatic than any other Bayesian model.

[Andrew Gelman]

We _do_ want to set the world ablaze. There’s lots of talk about black box VB but this would be the first truly black box implementation. This could be huge.

[Alp Kucukelbir]

> Depends what you want to do. These will drop constant terms and apply
> the Jacobian.

OK, I think that should do the trick. I'll ping again if confused. (Hopefully not.)

>> Finally: re being a good citizen, not setting the world ablaze, etc... what is it exactly that I should be doing as I implement variational bayes?
>
> I'm not sure what you're asking.

I just got worried based on the end of the meeting yesterday that I was committing some terrible crime as we speak. I didn't quite understand why you pointed me to stan/src/test/unit/model/util_test.cpp afterwards.

> Did you see the Developer Process
> wiki page on stan-dev/stan? Follow that. We're happy to clarify if
> it's unclear.

I found something on unit testing here: https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan

But I'm not sure how to do this for such a big thing like VB. I've written for instance a helper class called "stan::vb::latent_vars.hpp". Perhaps I should try to develop unit tests for that

Or perhaps the more honest question is: could someone ELI5 (explain like I'm five) what I should do right now, before writing any more VB code?

[Alp Kucukelbir]

I'm all in for that kind of arson! :)

As for ARD: I am certainly not married to that terminology. I was simply following Jan Drugowitsch's very nice tutorial-esque treatment of variational bayes for lin/log regression (http://arxiv.org/pdf/1310.5438.pdf).

[Bob Carpenter]

On Jul 25, 2014, at 3:32 PM, Alp Kucukelbir <akucu...@gmail.com> wrote:

>
>> Depends what you want to do. These will drop constant terms and apply
>> the Jacobian.
>
> OK, I think that should do the trick. I'll ping again if confused. (Hopefully not.)
>
>>> Finally: re being a good citizen, not setting the world ablaze, etc... what is it exactly that I should be doing as I implement variational bayes?
>>
>> I'm not sure what you're asking.
>
> I just got worried based on the end of the meeting yesterday that I was committing some terrible crime as we speak. I didn't quite understand why you pointed me to stan/src/test/unit/model/util_test.cpp afterwards.

Because that's where the test should go to show that
the two things you were trying to compute produced the
same answers. The stan/model/util.hpp file is where
the stan::model::gradient function is defined.

>
>> Did you see the Developer Process
>> wiki page on stan-dev/stan? Follow that. We're happy to clarify if
>> it's unclear.
>
> I found something on unit testing here: https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan
>
> But I'm not sure how to do this for such a big thing like VB. I've written for instance a helper class called "stan::vb::latent_vars.hpp". Perhaps I should try to develop unit tests for that

Yes. The point is that you should be developing small
pieces of functionality bottom-up with clean interfaces that
can be tested as you go. Those are the unit tests.

Then we'll put everything together and do functional tests
that the whole thing works end-to-end fitting models.

> Or perhaps the more honest question is: could someone ELI5 (explain like I'm five) what I should do right now, before writing any more VB code?

I'm not sure what you're asking. Just follow the developer
process for form, and building small, modular, testable units
for function.

- Bob

[Andrew Gelman]

Yes, to be clear we should use all terminology, we should call it hierarchical/multilevel modeling and also refer to the tutorial and say it’s also ARD.

[Alp Kucukelbir]

> Because that's where the test should go to show that
> the two things you were trying to compute produced the
> same answers. The stan/model/util.hpp file is where
> the stan::model::gradient function is defined.

Oh, I see! OK, so I need to create a folder stan/src/test/unit/vb and start writing my unit tests there.

> Yes. The point is that you should be developing small
> pieces of functionality bottom-up with clean interfaces that
> can be tested as you go. Those are the unit tests.

OK I will try writing and running some unit tests soon.

> Then we'll put everything together and do functional tests
> that the whole thing works end-to-end fitting models.

Very well. Sounds like a plan.

> I'm not sure what you're asking. Just follow the developer
> process for form, and building small, modular, testable units
> for function.

I will try my best. Sorry for theses adjustment pains... the learning curve for developing for Stan is non-trivial!

Cheers,
Alp

[Bob Carpenter]

On Jul 25, 2014, at 3:48 PM, Alp Kucukelbir <akucu...@gmail.com> wrote:

> ...

> Oh, I see! OK, so I need to create a folder stan/src/test/unit/vb and start writing my unit tests there.

There should be a test file for each source file with matching names.

> ...

> OK I will try writing and running some unit tests soon.

Just look at how the other tests are written --- it's pretty
easy with the googletest framework.

> ...

> I will try my best. Sorry for theses adjustment pains...
> the learning curve for developing for Stan is non-trivial!

Sorry there's not much doc on the code itself.

Our process is pretty typical in terms of how we're using
GitHub for branching, pull requests for code review, and
continuous integration testing.

- Bob