Possible Language Features

[Michael Betancourt]

Three questions about the validity of new features -- feel free to split up into separate threads when replying.

1. Prior/Likelihood Blocks

Why did we go with a single model block instead of separate prior/likelihood blocks? In order to make it easier to code models where the distinction between the two is vague? For philosophical reasons?

Technically it would be easy to add prior and likelihood blocks (while keeping a model block -- the final log_prob would just be model + prior + likelihood so that prior/likelhood are optional) which would make nested sampling, posterior predictive checks, etc, easier, no?

2. Non-differentiable Functions

I opened a feature request for exposing FFT to the language, and within a few minutes I realized how hard it would be due to the necessary gradients. But having an FFT for data preprocessing in the "transformed data" block would be extremely useful for many physics applications.

In terms of the language construction, is it difficult to limit functions to the "transformed data" block?

3. Jacobians

Looking through the transformation code, am I correct in assuming that

type_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y)

is called whenever transforming, unless the transformation is called within a probability density in which case

type_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y, T &lp)

is called? How does the language identify when a transformation is in the scope of a density?


In any case, right now we force the user to code in the Jacobians by hand but that seems extremely wasteful given that we already have them coded up. Can we expose

type_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y, T &lp)

into the language so that users can just call something like "log_jacobian_uniform(a, b)" to increment using the existing code? Refactoring all the transformations to isolate the Jacobians from the transformations would be a pain, but even calling "type_constrain" directly and accepting the overhead would probably be worth the ability to avoid coding the Jacobian by hand.

[Marcus Brubaker]

On Sun, Aug 18, 2013 at 2:09 PM, Michael Betancourt <betan...@gmail.com> wrote:

Three questions about the validity of new features -- feel free to split up into separate threads when replying.

1. Prior/Likelihood Blocks

Why did we go with a single model block instead of separate prior/likelihood blocks?  In order to make it easier to code models where the distinction between the two is vague?  For philosophical reasons?

Technically it would be easy to add prior and likelihood blocks (while keeping a model block -- the final log_prob would just be model + prior + likelihood so that prior/likelhood are optional) which would make nested sampling, posterior predictive checks, etc, easier, no?

I don't know why the decision was made originally, but in some models the division between likelihood and prior is rather arbitrary.  Also, on a practical level, there can be common terms or calculations between the prior and likelihood which can be shared with everything in one context.  That said, this is a useful distinction that I've been thinking about.  One possible solution I thought of would be to allow named sub-blocks in the model.  E.G., one could do:

model {

  likelihood {

    x ~ normal(mu,sigma);

  }

  prior {

    mu ~ normal(0,1);

    sigma ~ exponential(1);

  }

}

2. Non-differentiable Functions

I opened a feature request for exposing FFT to the language, and within a few minutes I realized how hard it would be due to the necessary gradients.  But having an FFT for data preprocessing in the "transformed data" block would be extremely useful for many physics applications.

In terms of the language construction, is it difficult to limit functions to the "transformed data" block?

What gradients in FFT are hard?  The FFT is just a linear transform, so its gradient is easy.  The hardest part of the FFT is that Stan has no way to deal with complex numbers right now.

 

3. Jacobians

Looking through the transformation code, am I correct in assuming that

type_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y)

is called whenever transforming, unless the transformation is called within a probability density in which case

type_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y, T &lp)

is called?  How does the language identify when a transformation is in the scope of a density?

The log_prob function has a template parameter which determines whether the volume correction is used.  Optimization doesn't use it, sampling does.  The version without lp is also called for I/O I believe.

 

In any case, right now we force the user to code in the Jacobians by hand but that seems extremely wasteful given that we already have them coded up.  Can we expose

type_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y, T &lp)

into the language so that users can just call something like "log_jacobian_uniform(a, b)" to increment using the existing code?  Refactoring all the transformations to isolate the Jacobians from the transformations would be a pain, but even calling "type_constrain" directly and accepting the overhead would probably be worth the ability to avoid coding the Jacobian by hand.

Seems reasonable and easy enough to do, but I'm not sure how useful it would be since if the transformation was correct wouldn't they have just used it directly?  Or is there some use case I'm missing?

Cheers,

Marcus

[Michael Betancourt]

> model {
> likelihood {
> x ~ normal(mu,sigma);
> }
> prior {
> mu ~ normal(0,1);
> sigma ~ exponential(1);
> }
> }

Yeah, any way so that underneath the hood you get three functions and gradients: prior, likelihood, and posterior.

> What gradients in FFT are hard? The FFT is just a linear transform, so its gradient is easy.

True -- I've been spending too much time with iterative ODE solvers lately.

> The hardest part of the FFT is that Stan has no way to deal with complex numbers right now.

Getting two vectors out and letting users do their own complex manipulations would be more than adequate.

> Seems reasonable and easy enough to do, but I'm not sure how useful it would be since if the transformation was correct wouldn't they have just used it directly? Or is there some use case I'm missing?

User-defined distributions using increment_logp where the Jacobians have to be added by hand.

[Marcus Brubaker]

On Sun, Aug 18, 2013 at 2:57 PM, Michael Betancourt <betan...@gmail.com> wrote:

> model {
>   likelihood {
>     x ~ normal(mu,sigma);
>   }
>   prior {
>     mu ~ normal(0,1);
>     sigma ~ exponential(1);
>   }
> }

Yeah, any way so that underneath the hood you get three functions and gradients: prior, likelihood, and posterior.

That could be harder in this construct because what if you had:

model {

  real tmp_x;

  likelihood {

    tmp_x <- 2 * x;

    x ~ something(tmp_x);

  }

  prior { 

   sigma ~ something_else(tmp_x);

  }

}

I was thinking more in terms of a single function which computed three values all at once, lp_prior, lp_likelihood and lp_posterior.  lp_posterior would be the sum of the other two plus the volume correction terms and anything that was done in a sub-block.

> The hardest part of the FFT is that Stan has no way to deal with complex numbers right now.

Getting two vectors out and letting users do their own complex manipulations would be more than adequate.

True, but then we'd need the ability to return more than one value from a function which isn't something that's done in the language right now.

 

> Seems reasonable and easy enough to do, but I'm not sure how useful it would be since if the transformation was correct wouldn't they have just used it directly?  Or is there some use case I'm missing?

User-defined distributions using increment_logp where the Jacobians have to be added by hand.

But if they defined the distribution on the unconstrained parameters then they don't need the Jacobian.  Alternatively, if they defined the distribution on the constrained parameters then why not just use a constrained definition?

Cheers,

Marcus

[Andrew Gelman]

I agree with Marcus.  I think that partitioning the factors of the model can be important, but this partitioning can be done in different ways.  A hierarchical model, for example, can be partitioned in more than one way.

For posterior predictive checks, the key point is that the user must decide what aspects of the model will persist and what aspects will be replicated.  This really has do be done as part of a graphical model structure, although I suspect there's a way to hack it in the current version of Stan.  But it goes deeper than prior/likelihood/posterior.

A

[Michael Betancourt]

> I was thinking more in terms of a single function which computed three values all at once, lp_prior, lp_likelihood and lp_posterior. lp_posterior would be the sum of the other two plus the volume correction terms and anything that was done in a sub-block.

Sure, but you could just as easily have three wrappers:

prior { return lp_prior + Jacobian correction }
likelihood { return lp_likelihood + Jacobian correction }
posterior { return lp_prior + lp_posterior + Jacobian correction }

> Getting two vectors out and letting users do their own complex manipulations would be more than adequate.
>
> True, but then we'd need the ability to return more than one value from a function which isn't something that's done in the language right now.

Or just return a 2xn matrix.

> But if they defined the distribution on the unconstrained parameters then they don't need the Jacobian. Alternatively, if they defined the distribution on the constrained parameters then why not just use a constrained definition?

There are lots of complex distributions in physics that are naturally defined on contained parameters. Yes, you can define them on unconstrained parameters but that is equivalent to computing and coding in the Jacobian. The idea is to make it much easier for the user by not requiring this additional step.

[Michael Betancourt]

> I agree with Marcus. I think that partitioning the factors of the model can be important, but this partitioning can be done in different ways. A hierarchical model, for example, can be partitioned in more than one way.

Agreed, which is why we don't want to force a prior/likelihood separation but offer it as an option for people implementing models where they have defined a separation.

[Andrew Gelman]

If so, I like the option of allowing people to give arbitrary names to different segments of their log-posterior. Of course the next step is to allow these to be a vector . . .

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

[Matt Hoffman]

I agree with Andrew on this. From the perspective of pure, ideal,
perfect, computation-free posterior inference, there's no distinction
between prior and likelihood, they're just independent terms in the
log-posterior. And we've been talking about a variety of situations
where we'd like to be able to partition the pieces of the
model—besides the PPCs and nested sampling applications mentioned
above, if we eventually want to be able to do stochastic inference
(e.g. contrastive divergence, stochastic EM, stochastic variational
inference) being able to partition the parameters (e.g. into point
estimable and nuisance parameters) would be really helpful.

Matt

[Michael Betancourt]

Agreed. My current wish list is

1) Higher-order autodiff
2) Parameter tagging
3) Prior/Likelihood Identification

Lots of fun stuff on the horizon.

[Bob Carpenter]

> On Aug 19, 2013, at 7:44 PM, Matt Hoffman <mdho...@cs.princeton.edu> wrote:
>
>> I agree with Andrew on this. From the perspective of pure, ideal,
>> perfect, computation-free posterior inference, there's no distinction
>> between prior and likelihood, they're just independent terms in the
>> log-posterior. And we've been talking about a variety of situations
>> where we'd like to be able to partition the pieces of the

>> model�besides the PPCs and nested sampling applications mentioned

>> above, if we eventually want to be able to do stochastic inference
>> (e.g. contrastive divergence, stochastic EM, stochastic variational
>> inference) being able to partition the parameters (e.g. into point
>> estimable and nuisance parameters) would be really helpful.

On 8/19/13 4:31 PM, Michael Betancourt wrote:
> Agreed. My current wish list is
>
> 1) Higher-order autodiff

This is our next big project after 2.0. We don't think
it'll be that hard other than the Windows compile issues.
It's already there for everything but the probability
functions. So if you write your model by hand, you can
use it now.

> 2) Parameter tagging

Agreed. This is most general. This should be pretty easy.

> 3) Prior/Likelihood Identification

Identification how? As in different functions in the model,
like if theta is the parameters, defining

log_prior(theta,...)

and

log_likelihood(theta,...);

so that

log_prob(theta,...) = log_prior(theta,...) + log_likelihood(theta,...).

I'm not sure we want every model to have these.

If we had statement tagging, we could identify the pieces in
the model, but I don't see how then to expose the different
functions. Maybe:

log_sub_prob(TAG,theta,...)

?

That would be general enough, I think, and probably also not too
difficult.

> Lots of fun stuff on the horizon.

Indeed!

(More below)

>> On Sun, Aug 18, 2013 at 12:43 PM, Andrew Gelman
>> <gel...@stat.columbia.edu> wrote:
>>> If so, I like the option of allowing people to give arbitrary names to different segments of their log-posterior. Of course the next step is to allow these to be a vector . . .
>>>
>>> On Aug 18, 2013, at 9:38 PM, Michael Betancourt wrote:
>>>
>>>>> I was thinking more in terms of a single function which computed three values all at once, lp_prior, lp_likelihood and lp_posterior. lp_posterior would be the sum of the other two plus the volume correction terms and anything that was done in a sub-block.
>>>>
>>>> Sure, but you could just as easily have three wrappers:
>>>>
>>>> prior { return lp_prior + Jacobian correction }
>>>> likelihood { return lp_likelihood + Jacobian correction }
>>>> posterior { return lp_prior + lp_posterior + Jacobian correction }
>>>>
>>>>> Getting two vectors out and letting users do their own complex manipulations would be more than adequate.
>>>>>
>>>>> True, but then we'd need the ability to return more than one value from a function which isn't something that's done in the language right now.
>>>>
>>>> Or just return a 2xn matrix.
>>>>
>>>>> But if they defined the distribution on the unconstrained parameters then they don't need the Jacobian. Alternatively, if they defined the distribution on the constrained parameters then why not just use a constrained definition?
>>>>
>>>> There are lots of complex distributions in physics that are naturally defined on contained parameters. Yes, you can define them on unconstrained parameters but that is equivalent to computing and coding in the Jacobian. The idea is to make it much easier for the user by not requiring this additional step.

I think we can do this. It's not on my near-term to-do list
given everything else we have going. We should actually be implementing
part of what we need with the ode solvers, which are going to need
Jacobians computed for arbitrary functions.

What we need is a way to specify it in Stan.

- Bob

[Michael Betancourt]

> Identification how? As in different functions in the model,
> like if theta is the parameters, defining
>
> log_prior(theta,...)
>
> and
>
> log_likelihood(theta,...);
>
> so that
>
> log_prob(theta,...) = log_prior(theta,...) + log_likelihood(theta,...).
>
> I'm not sure we want every model to have these.
>
> If we had statement tagging, we could identify the pieces in
> the model, but I don't see how then to expose the different
> functions. Maybe:
>
> log_sub_prob(TAG,theta,...)
>
> ?
>
> That would be general enough, I think, and probably also not too
> difficult.

That would be the most general. But having methods (even if they're wrappers)

prior { return lp_prior + Jacobian correction }
likelihood { return lp_likelihood + Jacobian correction }
posterior { return lp_prior + lp_posterior + Jacobian correction }

would be my ideal use case.

Again, posterior is not necessary equal to prior + likelihood. The idea would be to either have sub-blocks

model {
prior{
}
likelihood {
}
}

or three different blocks

model {}
prior {}
likelihood {}

with posterior = model + prior + likelihood so that the prior/likelihood terms are completely optional and would only be used if the user had decomposition in mind. Bonus points for being able to use the model/prior/likelihood functions in the generated quantities block.

My original motivation was that adding new blocks would just amount to duplicated the existing model parsing code.

[Bob Carpenter]

On 8/18/13 3:43 PM, Andrew Gelman wrote:
> If so, I like the option of allowing people to give arbitrary names to different segments of their log-posterior.

An issue we might have with this is how to calculate
any shared terms across these pieces (such as constants,
or a prior parameter that shows up in the prior and likelihood
functions).

- Bob

[Bob Carpenter]

On 8/18/13 2:37 PM, Marcus Brubaker wrote:

> On Sun, Aug 18, 2013 at 2:09 PM, Michael Betancourt <betan...@gmail.com <mailto:betan...@gmail.com>> wrote:
>
> Three questions about the validity of new features -- feel free to split up into separate threads when replying.
>
> 1. Prior/Likelihood Blocks
>
> Why did we go with a single model block instead of separate prior/likelihood blocks? In order to make it easier to
> code models where the distinction between the two is vague? For philosophical reasons?
>
> Technically it would be easy to add prior and likelihood blocks (while keeping a model block -- the final log_prob
> would just be model + prior + likelihood so that prior/likelhood are optional) which would make nested sampling,
> posterior predictive checks, etc, easier, no?
>
>
>
> I don't know why the decision was made originally, but in some models the division between likelihood and prior is
> rather arbitrary. Also, on a practical level, there can be common terms or calculations between the prior and
> likelihood which can be shared with everything in one context. That said, this is a useful distinction that I've been
> thinking about. One possible solution I thought of would be to allow named sub-blocks in the model. E.G., one could do:
>
> model {
> likelihood {
> x ~ normal(mu,sigma);
> }
> prior {
> mu ~ normal(0,1);
> sigma ~ exponential(1);
> }
> }

I'd still like a more general way to do this. Maybe something like
labeling statements:

@prior
alpha <- rep(0.0, 10);

@prior
theta_raw ~ foo(alpha);

@likelihood
theta <- bar(theta_raw);

@likelihood
y ~ mulitnomial(theta);

As Michael's pointed out, terms may be shared, so we could do something
like:

@prior @likelihood
z <- baz(theta);

> 2. Non-differentiable Functions
>
> I opened a feature request for exposing FFT to the language, and within a few minutes I realized how hard it would
> be due to the necessary gradients. But having an FFT for data preprocessing in the "transformed data" block would
> be extremely useful for many physics applications.
>
> In terms of the language construction, is it difficult to limit functions to the "transformed data" block?
>
>
>
> What gradients in FFT are hard? The FFT is just a linear transform, so its gradient is easy. The hardest part of the
> FFT is that Stan has no way to deal with complex numbers right now.

If someone can write out the derivatives, it should be easy to
code them up by hand.

But we don't need the complexes if the inputs are all real, do we?

> 3. Jacobians
>
> Looking through the transformation code, am I correct in assuming that
>
> type_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y)
>
> is called whenever transforming, unless the transformation is called within a probability density in which case
>
> type_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y, T &lp)
>
> is called? How does the language identify when a transformation is in the scope of a density?



> The log_prob function has a template parameter which determines whether the volume correction is used. Optimization
> doesn't use it, sampling does. The version without lp is also called for I/O I believe.

No, you still need lp for I/O because we write out lp__.

> In any case, right now we force the user to code in the Jacobians by hand but that seems extremely wasteful given
> that we already have them coded up. Can we expose
>
> type_constrain(const Eigen::Matrix<T,Eigen::Dynamic,1>& y, T &lp)

Absolutely. It's on the to-do list.

It's no longer trivial because we've deprecated lp__. So these
somehow need to get ahold of lp__ directly in the same way as
the sampling statements do. And the same way the RNG need an
implicit RNG. Shouldn't be too tricky.

> into the language so that users can just call something like "log_jacobian_uniform(a, b)" to increment using the
> existing code? Refactoring all the transformations to isolate the Jacobians from the transformations would be a
> pain, but even calling "type_constrain" directly and accepting the overhead would probably be worth the ability to
> avoid coding the Jacobian by hand.
>
>
>
> Seems reasonable and easy enough to do, but I'm not sure how useful it would be since if the transformation was correct
> wouldn't they have just used it directly? Or is there some use case I'm missing?

That's why it wasn't a high priority to add. That, and because
it's another sub-case to doc that might confuse users.

- Bob

[Michael Betancourt]

Wouldn't we assume that all parameters are accessible in each block?

As for constants and temps, why wouldn't we just treat the blocks as separate scopes?

[Bob Carpenter]

On 8/20/13 4:31 PM, Michael Betancourt wrote:
> Wouldn't we assume that all parameters are accessible in each block?

That would make sense. And all the data and transformed data.

> As for constants and temps, why wouldn't we just treat the blocks as separate scopes?

That's simplest, but it prevents sharing local variables across
blocks. And it forces all the transformed parameters to be computed
for each block if you want to also include those.

- Bob

[Michael Betancourt]

>> What gradients in FFT are hard? The FFT is just a linear transform, so its gradient is easy. The hardest part of the
>> FFT is that Stan has no way to deal with complex numbers right now.
>
> If someone can write out the derivatives, it should be easy to
> code them up by hand.
>
> But we don't need the complexes if the inputs are all real, do we?

You always get complex outputs -- equivalently frequency and phase info.

[Michael Betancourt]

> Wouldn't we assume that all parameters are accessible in each block?
>
> That would make sense. And all the data and transformed data.
>
>> As for constants and temps, why wouldn't we just treat the blocks as separate scopes?
>
> That's simplest, but it prevents sharing local variables across
> blocks. And it forces all the transformed parameters to be computed
> for each block if you want to also include those.

I guess what I had imagined is that we don't literally have

log_post() { return log_prior() + log_like(); }

but rather have the parser dump the log_prior and log_like code into log_post (along with any code especially specified in the model block), so that transformed parameters would not have to be computed multiple times.

I also think not being able to share local variables is a fair sacrifice for the ability to specify prior/likelihood, but I may be in the minority.