forward-looking transformations

[Ben Goodrich]

So, we have a several examples of Matt-trickish things where the optimal thing to do would be to put a prior on primitive parameters and transform them to imply a particular prior on substantive parameters. Currently, we are forcing users to do all these tricks in the .stan file, which can be tedious and / or error prone. We could write convenience functions to do the transformations (and maybe optimize the derivatives internally), but it would still be up to the user to put the appropriate prior on the primitive parameters, which would still be tedious and error-prone.

Is there some way we can design things so that particular transformation functions get to change the log-posterior? So, something like

parameters {
  vector[K] z;
}
model {
  vector[K] s;
  s <- normalize2uniform(z);
}

would be equivalent to

parameters {
  vector[K] z;
}
model {
  vector[K] s;
  s <- z / sqrt(dot_self(z))
  z ~ normal(0.0, 1.0);
}

I guess either Stan would need to know that it needs to pass the log-posterior by value to functions like normalize2uniform or we would have to rearrange the scope so that such functions could modify lp__ .

Ben

[Bob Carpenter]

On 11/15/13, 5:35 PM, Ben Goodrich wrote:
> So, we have a several examples of Matt-trickish things where the optimal thing to do would be to put a prior on
> primitive parameters and transform them to imply a particular prior on substantive parameters. Currently, we are forcing
> users to do all these tricks in the .stan file, which can be tedious and / or error prone.

Agreed.

> We could write convenience
> functions to do the transformations (and maybe optimize the derivatives internally), but it would still be up to the
> user to put the appropriate prior on the primitive parameters, which would still be tedious and error-prone.

But in this case, those would be standardized and
could get matching functions or it could all be rolled
into a single function.

> Is there some way we can design things so that particular transformation functions get to change the log-posterior? So,

Yes, it's possible. It's basically how the
log probability functions. The parser recognizes
a function ending in "_log" (a very bad choice in
retrospect) and adds an argument for lp__, which
is passed by non-constant reference.

Scope of where these functions are allowed is restricted
to the blocks where there is a log probability to
modify --- transformed parameters and model.

> something like
>
> |
> parameters {
> vector[K] z;
> }
> model {
> vector[K] s;
> s <- normalize2uniform(z);
> }
> |
>
> would be equivalent to
>
> |
> parameters {
> vector[K] z;
> }
> model {
> vector[K] s;
> s <- z / sqrt(dot_self(z))
> z ~ normal(0.0, 1.0);
> }

I'm afraid I can't quite follow the intent of this
transform. What's going to happen with the s?

But whatever it's intent, we could code the
first thing so that it behaves like the second in just
the way you suggest.

I was thinking we'd transform something like:

parameters {
real mu;
real sigma;
real y;
}
model {
y ~ normal(mu,sigma);
}

to be equivalent to this:

parameters {
real mu;
real sigma;
real y_unit;
}
transformed parameters {
real y;
y <- mu + sigma * y_unit;
}
model {
y_unit ~ normal(0,1);
}

But that I'm not so sure how to do because it's not
just adding a function, and it drops an implicit y_raw
value on users. And I believe there will be ordering
issues that are complicated if there are multiple instances.

> I guess either Stan would need to know that it needs to pass the log-posterior by value to functions like
> normalize2uniform or we would have to rearrange the scope so that such functions could modify lp__ .

I think you meant by non-constant reference rather than by
value. Most arguments are passed by constant reference. All
we have to do is drop the "const" in the definition and it can
modify its inputs. So far, I've tried to stay away from doing
that explicitly because I think reference vs. value is a complicated
enough concept that it'll confuse users. Passing by const
reference gives us the efficiency of pass by ref with the
transparency of pass-by-value.

- Bob

[Michael Betancourt]

I am less hopeful on automated transformations. The problem is that non-centered
parameterizations are not always the best approach (so you'll have to enable some
ability to turn any automated transformation on and off) and in many cases they
don't even make sense (anything beyond a one-level hierarchy). Personally I don't
think locking users into building one-layer multilevel models is worth the possible
convenience.

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

[Bob Carpenter]

That doesn't prohibit building one level at a time
in a special function or distribution. Anything we
can do to more encapsulate the process would probably
help.

- Bob

[Michael Betancourt]

The best you'd be able to do is encapsulate a one-level model as a model template
or the like. If one can figure out the notation for writing a "fast one-way" model then sure,
but this can't be used as a building block to bigger models because
the default non-centered parameterizations are no longer valid when you have multiple
layers (in particular, you can only "non-center" one layer in a hierarchy, and which
to non-center can be extremely difficult to determine).

[Ben Goodrich]

On Friday, November 15, 2013 6:29:51 PM UTC-5, Bob Carpenter wrote:

> Is there some way we can design things so that particular transformation functions get to change the log-posterior? So,

Yes, it's possible.  It's basically how the
log probability functions.  The parser recognizes
a function ending in "_log" (a very bad choice in
retrospect) and adds an argument for lp__, which
is passed by non-constant reference.

Okay, so these convenience functions in C++ just need to have the magic _log suffix? That is a bit hackish but will work fine. Whether we want to ever do automatic transformations is a separate debate, but having convenient functions for users to call seems unambiguously good. I'll try one on my onion branch.

Ben

[Bob Carpenter]

Agreed. There may already be an issue for adding the
transforms as functions --- it's definitely on the to-do list.

We can come up with another kind of convention or instead
make a list and keep it in a set-like data structure.

We're also going to have to think about better error reporting
for functions like sqrt() that are in the lib. If we want to pass
in code location for debugging to all of our calls, then we could
also just always pass in lp__ as the last arg.

- Bob

[Ben Goodrich]

On Sunday, November 17, 2013 11:29:20 AM UTC-5, Bob Carpenter wrote:

On 11/15/13, 8:48 PM, Ben Goodrich wrote:
> On Friday, November 15, 2013 6:29:51 PM UTC-5, Bob Carpenter wrote:
>
>      > Is there some way we can design things so that particular transformation functions get to change the
>     log-posterior? So,
>
>     Yes, it's possible.  It's basically how the
>     log probability functions.  The parser recognizes
>     a function ending in "_log" (a very bad choice in
>     retrospect) and adds an argument for lp__, which
>     is passed by non-constant reference.
>
>
> Okay, so these convenience functions in C++ just need to have the magic _log suffix? That is a bit hackish but will work
> fine. Whether we want to ever do automatic transformations is a separate debate, but having convenient functions for
> users to call seems unambiguously good. I'll try one on my onion branch.

Agreed.  There may already be an issue for adding the
transforms as functions --- it's definitely on the to-do list.

I think we may have created an issue for implementing commonly-used Jacobians, which would be another case where the function needs to be able to change lp__ . In the cases I was referring to, the function needs to be able to change lp__ but there is no explicit Jacobian to calculate in the code.
 

We can come up with another kind of convention or instead
make a list and keep it in a set-like data structure.

That might be better.
 

We're also going to have to think about better error reporting
for functions like sqrt() that are in the lib.  If we want to pass
in code location for debugging to all of our calls, then we could
also just always pass in lp__ as the last arg.

Always pass it in but for most functions make it a const reference? I guess that could work but it seems like there should be a better way.

Ben

[Bob Carpenter]

On 11/17/13, 4:54 PM, Ben Goodrich wrote:
> On Sunday, November 17, 2013 11:29:20 AM UTC-5, Bob Carpenter wrote:
>
> On 11/15/13, 8:48 PM, Ben Goodrich wrote:
> > On Friday, November 15, 2013 6:29:51 PM UTC-5, Bob Carpenter wrote:
> >
> > > Is there some way we can design things so that particular transformation functions get to change the
> > log-posterior? So,
> >
> > Yes, it's possible. It's basically how the
> > log probability functions. The parser recognizes
> > a function ending in "_log" (a very bad choice in
> > retrospect) and adds an argument for lp__, which
> > is passed by non-constant reference.
> >
> >
> > Okay, so these convenience functions in C++ just need to have the magic _log suffix? That is a bit hackish but
> will work
> > fine. Whether we want to ever do automatic transformations is a separate debate, but having convenient functions for
> > users to call seems unambiguously good. I'll try one on my onion branch.
>
> Agreed. There may already be an issue for adding the
> transforms as functions --- it's definitely on the to-do list.
>
>
> I think we may have created an issue for implementing commonly-used Jacobians, which would be another case where the
> function needs to be able to change lp__ . In the cases I was referring to, the function needs to be able to change lp__
> but there is no explicit Jacobian to calculate in the code.

OK.

> We can come up with another kind of convention or instead
> make a list and keep it in a set-like data structure.
>
>
> That might be better.

> We're also going to have to think about better error reporting
> for functions like sqrt() that are in the lib. If we want to pass
> in code location for debugging to all of our calls, then we could
> also just always pass in lp__ as the last arg.
>
>
> Always pass it in but for most functions make it a const reference? I guess that could work but it seems like there
> should be a better way.

We can do what Boost static config does and actually define
the functions with no arguments, just the types (which you
need to resolve which function is called). In a typical
templated function for us, that'd look like this:

template <typename T1, ..., typename TN>
inline
typename promote_args<T1,...,TN>::type
foo(const typename promote_args<T1,...,TN>::type& /* ignore lp arg */,
int src_line,

... other args ...) {

... do something ...

... if something goes wrong, report source position ...
}

We could make the source_line a size_t type and we could also
include a source column.

This is a pretty big change, but we can tackle it function by function
if we get the right conventions into place.

Both the testing and the passing of source and the deterrent to
inlining by making the functions bigger and more conditional will
potentially cause a noticeable slowdown.

So if we still want to have a way of turning the error checking
off, that's much trickier. Then we have to dummy out the source_line
arg with an ifdef along with the tests.

- Bob

[Marcus Brubaker]

What if we just provide the transform and Jacobian determinant functions as normal functions?  E.G.,  allow users to write

parameters {

  vector[N] z;

}

transformed parameters {

  vector[N] positive_z;

  positive_z <- exp_transform(z);

  increment_log_prob(exp_transform_lpcorrection(z));

}

This way exp_transform() and exp_transform_lpcorrection() are just normal functions.

Alternatively, maybe we use a different operator for this?  Maybe =  E.G.,

  

transformed parameters {

  vector[N] positive_z;

  positive_z = exp_transform(z);

}

Which then implicitly includes the increment_log_prob(exp_transform_lpcorrection(z)) line?

The nice side-effect of this notation (or any other similar one) is, when we can do higher order autodiff or symbolic differentiation, we could use this for more general transforms and then automatically compute the Jacobian of the transform.

What if we do what we've done with some of the transformation code: add try{}catch{} blocks to the generated code which catches exceptions, prepends a message with the the source line number and rethrows it.  Or, if we wanted to avoid many individual try{}catch{} blocks we could do something like:

int source_line_number = 0;

try {

  source_line_number = 1;

  ... do line 1 stuff ...

  source_line_number = 4;

  ... do line 4 stuff ...

}

catch () {

  ... report catching exception that was generated on line number source_line_number....

}

This might be a bit cleaner.  In either case, if we wanted to be really fancy, it would be nice if we could even print the model code along with it's line number.

I'd looked into doing something like this in the past, and the issue has always been getting the line numbers out of the parser and keeping them around when generating the code.  I'd never been able to penetrate the parser code enough to figure out how to do it.

Cheers,

Marcus

[Michael Betancourt]

I always thought we were doing the Jacobians with function-like calls as Marcus described below,
was I misinformed?

[Marcus Brubaker]

There are a bunch of constrain/unconstrain functions which implement the transforms and increment the log_prob but these aren't exposed to users and aren't really in a form that could be easily exposed to the modelling language.

[Bob Carpenter]

They're in prob/transform.hpp and should be broken
out into individual files, one per transform function.

I don't think they'd be that hard to add if we had
a way to pass in the lp__ variable from the modeling
language.

An orthogonal high-priority to-do item is to write
custom derivatives for the transforms. And to
vectorize them.

- Bob

On 11/18/13, 12:33 PM, Marcus Brubaker wrote:
> There are a bunch of constrain/unconstrain functions which implement the transforms and increment the log_prob but these
> aren't exposed to users and aren't really in a form that could be easily exposed to the modelling language.
>
>
> On Mon, Nov 18, 2013 at 12:27 PM, Michael Betancourt <betan...@gmail.com <mailto:betan...@gmail.com>> wrote:
>
> I always thought we were doing the Jacobians with function-like calls as Marcus described below,
> was I misinformed?
>

> On Nov 18, 2013, at 6:02 PM, Marcus Brubaker <marcus.a...@gmail.com <mailto:marcus.a...@gmail.com>> 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>.

[Daniel Lee]

+1 for breaking apart the prob/transform.hpp, writing custom derivatives, and vectorizing.

On Mon, Nov 18, 2013 at 1:03 PM, Bob Carpenter <ca...@alias-i.com> wrote:

They're in prob/transform.hpp and should be broken
out into individual files, one per transform function.

I don't think they'd be that hard to add if we had
a way to pass in the lp__ variable from the modeling
language.

An orthogonal high-priority to-do item is to write
custom derivatives for the transforms.  And to
vectorize them.

- Bob

On 11/18/13, 12:33 PM, Marcus Brubaker wrote:

There are a bunch of constrain/unconstrain functions which implement the transforms and increment the log_prob but these
aren't exposed to users and aren't really in a form that could be easily exposed to the modelling language.

On Mon, Nov 18, 2013 at 12:27 PM, Michael Betancourt <betan...@gmail.com <mailto:betan...@gmail.com>> wrote:

    I always thought we were doing the Jacobians with function-like calls as Marcus described below,
    was I misinformed?

    stan-dev+unsubscribe@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+unsubscribe@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+unsubscribe@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+unsubscribe@googlegroups.com.