Initialise particular parameters

[John Reid]

Hi, I have a model in which initialising a particular parameter vector helps it to find better posterior modes. However if I initialise this vector through the "init" argument to the Rstan stan() method, I think I have to initialise all the parameters. Is there anyway to avoid this? If not, could you add this as a feature? That is, if the user-supplied init function does not return a list element for a parameter then Stan would use "random" as the default.
Thanks,
John.

[Bob Carpenter]

Correct, it's either all or nothing for initialization.
Allowing partial inits has been on our radar, but we've
never gotten around to it.

I added a feature request to our issue tracker:

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

The random inits are by default [-2,2] on the constrained
scale and the transforms to constrained params are all
defined in a later chapter in the manual. It the transforms
were exposed in the interfaces, this would be easier, at
least in RStan and PyStan.

- Bob

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

[John Reid]

It is good to know this might turn up in a later Stan.

On a related note I was wondering if Stan is aware of the generative nature of some models. Most models I work with are described in terms of conditional distributions that have a natural order from hyperparameters through latent variables to data. I have had success initialising Stan using draws from the prior when the default "random" method didn't do so well. Are there any plans to incorporate some sort of initialisation like this? This would save users recoding their priors in init() functions (which can be error-prone).

Thanks,
John.

[Bob Carpenter]

On Oct 7, 2014, at 2:28 PM, John Reid <johnba...@gmail.com> wrote:

...

> On a related note I was wondering if Stan is aware of the generative nature of some models. Most models I work with are described in terms of conditional distributions that have a natural order from hyperparameters through latent variables to data. I have had success initialising Stan using draws from the prior when the default "random" method didn't do so well. Are there any plans to incorporate some sort of initialisation like this? This would save users recoding their priors in init() functions (which can be error-prone).

Short answers: No and no.

This has been something Andrew's been asking for, too.
But as is, Stan has no knowledge of any generative structure
in the model, so until we add a layer on top to deal with
directed graphical models or some more restricted class, we
don't have any hooks on which to hang such an implementation.

But I'd imagine in cases where what you're suggesting would
work, Stan would also be able to just converge.

Do you have non heavily multimodal examples where initialization
makes a big difference? If so, we'd love to collect them as
test cases.

My understanding is that JAGS just initializes from the prior
means and that BUGS does what you suggest, but I could be
misremembering or confusing the two systems.

- Bob

[John Reid]

On 07/10/14 22:17, Bob Carpenter wrote:
> On Oct 7, 2014, at 2:28 PM, John Reid <johnba...@gmail.com> wrote:
>
> ...
>
>> On a related note I was wondering if Stan is aware of the generative nature of some models. Most models I work with are described in terms of conditional distributions that have a natural order from hyperparameters through latent variables to data. I have had success initialising Stan using draws from the prior when the default "random" method didn't do so well. Are there any plans to incorporate some sort of initialisation like this? This would save users recoding their priors in init() functions (which can be error-prone).
> Short answers: No and no.
>
> This has been something Andrew's been asking for, too.
> But as is, Stan has no knowledge of any generative structure
> in the model, so until we add a layer on top to deal with
> directed graphical models or some more restricted class, we
> don't have any hooks on which to hang such an implementation.

I thought that would be the case.

There might be a half-way house though. If I code my model in Stan and
then write the initialisation function in R, I often have to transform
between different parameterisations. Perhaps there could be some sort of
"parameter initialisation" block in Stan models so this wouldn't be
necessary. Sharing user-defined functions between this block and the
model could be useful. Also if custom initialisation is so useful for
good inference it would be nice to include the initialisations used with
the definition of the model.

>
> But I'd imagine in cases where what you're suggesting would
> work, Stan would also be able to just converge.
>
> Do you have non heavily multimodal examples where initialization
> makes a big difference? If so, we'd love to collect them as
> test cases.

You've hit the nail on the head, my examples have multi-modal
posteriors. I believe initialising from the prior allows the sampler to
have a much greater chance of visiting more of the modes. I'm trying to
learn a latent temporal variable that explains the data well. Each datum
has a latent pseudotime and the model's parameters are allowed to
smoothly vary over this pseudotime. The initialisation of the data
pseudotime parameters seems very important as Stan has trouble
reordering the data during sampling. I guess that's not the sort of
example you're looking for.

Regards,
John.

[Daniel Lee]

On Wed, Oct 8, 2014 at 4:26 AM, John Reid <johnba...@gmail.com> wrote:

On 07/10/14 22:17, Bob Carpenter wrote:
> On Oct 7, 2014, at 2:28 PM, John Reid <johnba...@gmail.com> wrote:
>
> ...
>
>> On a related note I was wondering if Stan is aware of the generative nature of some models. Most models I work with are described in terms of conditional distributions that have a natural order from hyperparameters through latent variables to data. I have had success initialising Stan using draws from the prior when the default "random" method didn't do so well. Are there any plans to incorporate some sort of initialisation like this? This would save users recoding their priors in init() functions (which can be error-prone).
> Short answers: No and no.
>
> This has been something Andrew's been asking for, too.
> But as is, Stan has no knowledge of any generative structure
> in the model, so until we add a layer on top to deal with
> directed graphical models or some more restricted class, we
> don't have any hooks on which to hang such an implementation.
I thought that would be the case.

There might be a half-way house though. If I code my model in Stan and
then write the initialisation function in R, I often have to transform
between different parameterisations. Perhaps there could be some sort of
"parameter initialisation" block in Stan models so this wouldn't be
necessary. Sharing user-defined functions between this block and the
model could be useful. Also if custom initialisation is so useful for
good inference it would be nice to include the initialisations used with
the definition of the model.

Hi John, just curious: have you found the diffuse starting points to be a problem in your models? Stan initializes all the unconstrained parameters i.i.d. uniform between -2 and 2 inclusive. When initializing from R, it should be an initialization on the constrained space. That is, if a parameter is declared as `real<lower=0,upper=1> theta;`, then it should be initialized between (0, 1).

What transformations are you doing between different parameterizations? Hopefully you're not trying to go from the constrained to the unconstrained space from R outside of Stan.

>
> But I'd imagine in cases where what you're suggesting would
> work, Stan would also be able to just converge.
>
> Do you have non heavily multimodal examples where initialization
> makes a big difference?  If so, we'd love to collect them as
> test cases.
You've hit the nail on the head, my examples have multi-modal
posteriors. I believe initialising from the prior allows the sampler to
have a much greater chance of visiting more of the modes. I'm trying to
learn a latent temporal variable that explains the data well. Each datum
has a latent pseudotime and the model's parameters are allowed to
smoothly vary over this pseudotime. The initialisation of the data
pseudotime parameters seems very important as Stan has trouble
reordering the data during sampling. I guess that's not the sort of
example you're looking for.

It is an example of what we're looking for. Without knowing the details of your model: reordering is tough, multi-modality is tough. If specific initial values are putting Stan into a single mode, that's interesting to know and perhaps one of the future generations of Stan will be able to solve the problem without a hitch (I'm thinking Riemmanian manifold HMC). Of course, if it's really multi-modal and you're not trying to do a full Bayesian inference, that opens up a lot of doors that won't get you the right expectation, but might give you something from a high posterior mode.

 

Regards,

[Bob Carpenter]

On Oct 8, 2014, at 4:26 AM, John Reid <johnba...@gmail.com> wrote:

> On 07/10/14 22:17, Bob Carpenter wrote:

>> On Oct 7, 2014, at 2:28 PM, John Reid <johnba...@gmail.com> wrote:

> If I code my model in Stan and
> then write the initialisation function in R, I often have to transform
> between different parameterisations.

Different parameterizations of the model in Stan?

The inits are always on the constrained scale.

> Perhaps there could be some sort of
> "parameter initialisation" block in Stan models so this wouldn't be
> necessary.

This is something that we had in the works before our adaptation
got a lot more robust. We may still do it at some point.

> Sharing user-defined functions between this block and the
> model could be useful.

Agreed.

We plan to add the transform functions as built-ins soon along with all of
our internal error checking.

> Also if custom initialisation is so useful for
> good inference

If custom initialization is really necessary, it may signal
a problem with inference in the model as Daniel pointed out.

> it would be nice to include the initialisations used with
> the definition of the model.

Definitely.

>> But I'd imagine in cases where what you're suggesting would
>> work, Stan would also be able to just converge.
>>
>> Do you have non heavily multimodal examples where initialization
>> makes a big difference? If so, we'd love to collect them as
>> test cases.

> You've hit the nail on the head, my examples have multi-modal
> posteriors. I believe initialising from the prior allows the sampler to
> have a much greater chance of visiting more of the modes.

That will depend on whether unconstraining uniform(-2,2) or
the prior is more dispersed and other factors about the posterior
curvature.

In any case, you may visit more, but it's going to be impossible
to do proper Bayesian inference in highly multimodal posteriors.

We'll be adding Michael Betancourt's thermodynamic sampling (arXiv-ed),
which will probably help with a lot of these issues.

We may also may wind up doing some kind of tempered optimization,
which may be more relevant for inference under multiple modes.

> I'm trying to
> learn a latent temporal variable that explains the data well. Each datum
> has a latent pseudotime and the model's parameters are allowed to
> smoothly vary over this pseudotime. The initialisation of the data
> pseudotime parameters seems very important as Stan has trouble
> reordering the data during sampling. I guess that's not the sort of
> example you're looking for.

These are the sorts of examples we're looking for.
Unfortunately, anything with modes separated by low density
troughs is going to be pretty much impossible to mix in
practice. Like the mixture models in the manual: K-means
clustering, HMMs, and LDA. And even if you could mix, the
combinatorics will prevent proper posterior inference unless
one mode is so dominant the others can be ignored.

In this particular case, if the ordering is known, we have
a constrained parameter for ordered input.

- Bob

[Andrew Gelman]

Really?  I thought we’d swiched this to be uniform from -0.5 to 0.5 as default, as we were having problems with some of the more extreme starting points.

Also, in response to the earlier question, it would not be possible for Stan to automatically draw parameters from the prior (indeed “the prior” is not so clearly defined, as all we are specifying is a joint density function) but it would be possible to do this “manually” by doing such simulations in R and then specifying these as starting values in the Stan call.

A

[Andrew Gelman]

On Oct 8, 2014, at 1:28 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
>
> In any case, you may visit more, but it's going to be impossible
> to do proper Bayesian inference in highly multimodal posteriors.

“Impossible” is a bit strong. Let’s just say that default Nuts (or Metropolis, or Gibbs, or HMC) will have problems with multimodal posteriors.

> We'll be adding Michael Betancourt's thermodynamic sampling (arXiv-ed),
> which will probably help with a lot of these issues.

That’s “adiabatic sampling” now!

> Unfortunately, anything with modes separated by low density
> troughs is going to be pretty much impossible to mix in
> practice.

Again, just to clarify: “pretty much impossible to mix” using certain algorithms, maybe not with others.

[Bob Carpenter]

On Oct 8, 2014, at 4:14 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

> Really? I thought we’d swiched this to be uniform from -0.5 to 0.5 as default, as we were having problems with some of the more extreme starting points.

I thought so, too, but we never did.

The problems with extreme starting points went away with
better adaptation and higher acceptance rates, so we left it
at [-2,2].

It's easy to specify with init = 0.5 in RStan. The choice of 2
remains rather arbitrary. It covers a lot of logit-scaled parameters
and runs from roughly 0.1 to 10.0 on positive constrained parameters,
and we thought that looked suitably diffuse for a default.

> Also, in response to the earlier question, it would not be possible for Stan to automatically draw parameters from the prior (indeed “the prior” is not so clearly defined, as all we are specifying is a joint density function) but it would be possible to do this “manually” by doing such simulations in R and then specifying these as starting values in the Stan call.

Once we add the prior/likelihood (aka base/delta split), this
should be easy for models with proper priors.

Of course, it'll introduce all the confusion about prior/likelihood
we've been trying to avoid, but that's inevitable with the above split.

- Bob

[Bob Carpenter]

On Oct 8, 2014, at 4:16 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

>
> On Oct 8, 2014, at 1:28 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>>
>>
>> In any case, you may visit more, but it's going to be impossible
>> to do proper Bayesian inference in highly multimodal posteriors.
>
> "Impossible" is a bit strong.

When I say "highly multimodal", I'm talking about cases like LDA
or K-means clustering where the posterior multimodality is combinatorially
intractable (i.e., the number of modes grow exponentially and thus overwhelm any
hardware-based solution to explicitly sample from them).

> Let's just say that default Nuts (or Metropolis, or Gibbs, or HMC) will have problems with multimodal posteriors.

If someone could solve the combinatorial intractability problem,
we have bigger fish to fry than Bayesian sampling!

>> We'll be adding Michael Betancourt's thermodynamic sampling (arXiv-ed),
>> which will probably help with a lot of these issues.
>
> That's "adiabatic sampling" now!

So it is. Michael even changed the title of the arXiv paper:

http://arxiv.org/abs/1405.3489

>> Unfortunately, anything with modes separated by low density
>> troughs is going to be pretty much impossible to mix in
>> practice.
>
> Again, just to clarify: "pretty much impossible to mix" using certain algorithms, maybe not with others.

Right. Moving between a small number of modes is not intractable in
the same way that sampling from all of them is.

- Bob

[John Reid]

On 08/10/14 14:44, Daniel Lee wrote:
> Hi John, just curious: have you found the diffuse starting points to
> be a problem in your models? Stan initializes all the unconstrained
> parameters i.i.d. uniform between -2 and 2 inclusive. When
> initializing from R, it should be an initialization on the constrained
> space. That is, if a parameter is declared as `real<lower=0,upper=1>
> theta;`, then it should be initialized between (0, 1).

What I have found is that if I sample from the prior then I find more
and better modes in the posterior of my model. I understand inference in
multimodal models is difficult and I would be very happy just to find
the best mode rather than samples from the posterior.

>
> What transformations are you doing between different
> parameterizations? Hopefully you're not trying to go from the
> constrained to the unconstrained space from R outside of Stan.

I wasn't clear enough: I meant different parameterisations of a
distribution, for instance R might parameterise a distribution using a
transformation of the parameters that Stan uses. No I'm not trying to
initialise in the unconstrained space.

Regards,
John.

[John Reid]

On 08/10/14 18:28, Bob Carpenter wrote:
> On Oct 8, 2014, at 4:26 AM, John Reid <johnba...@gmail.com> wrote:
>
>> On 07/10/14 22:17, Bob Carpenter wrote:
>>> On Oct 7, 2014, at 2:28 PM, John Reid <johnba...@gmail.com> wrote:
>> If I code my model in Stan and
>> then write the initialisation function in R, I often have to transform
>> between different parameterisations.
> Different parameterizations of the model in Stan?
>
> The inits are always on the constrained scale.

Yes I understand that, I meant different parameterisations of a
distribution between Stan and R. E.g. the negative binomial has several
alternative parameterisations and switching between them in R and Stan
is a bit error prone. I would prefer to do almost everything model
related in Stan if I can.

>
>> Perhaps there could be some sort of
>> "parameter initialisation" block in Stan models so this wouldn't be
>> necessary.
> This is something that we had in the works before our adaptation
> got a lot more robust. We may still do it at some point.
>
>> Sharing user-defined functions between this block and the
>> model could be useful.
>
> Agreed.
>
> We plan to add the transform functions as built-ins soon along with all of
> our internal error checking.
>
>
>> Also if custom initialisation is so useful for
>> good inference
>
> If custom initialization is really necessary, it may signal
> a problem with inference in the model as Daniel pointed out.

Sure but as I understand it if that is my model I can either try custom
initialisation, reparameterise the model or try another inference
engine/method. Custom initialisation has given what I want so far.
....

>
>> I'm trying to
>> learn a latent temporal variable that explains the data well. Each datum
>> has a latent pseudotime and the model's parameters are allowed to
>> smoothly vary over this pseudotime. The initialisation of the data
>> pseudotime parameters seems very important as Stan has trouble
>> reordering the data during sampling. I guess that's not the sort of
>> example you're looking for.
> These are the sorts of examples we're looking for.
> Unfortunately, anything with modes separated by low density
> troughs is going to be pretty much impossible to mix in
> practice. Like the mixture models in the manual: K-means
> clustering, HMMs, and LDA. And even if you could mix, the
> combinatorics will prevent proper posterior inference unless
> one mode is so dominant the others can be ignored.

I'm happy just to find the best (or a good) mode. Full posterior
inference would be great but not necessary.

>
> In this particular case, if the ordering is known, we have
> a constrained parameter for ordered input.

The ordering is what I'm trying to infer. Again I wasn't quite clear, I
meant ordering in terms of the pseudotime. I'll see if I can tidy up the
model and post it here.

Thanks,
John.

[John Reid]

On 08/10/14 21:42, Bob Carpenter wrote:
>
> We'll be adding Michael Betancourt's thermodynamic sampling (arXiv-ed),
> which will probably help with a lot of these issues.
>> That's "adiabatic sampling" now!
>
> So it is. Michael even changed the title of the arXiv paper:
>
> http://arxiv.org/abs/1405.3489
>

Great, that looks very interesting.

Regards,
John.

[John Reid]

On Thursday, 9 October 2014 10:32:49 UTC+1, John Reid wrote:

On 08/10/14 18:28, Bob Carpenter wrote:

> On Oct 8, 2014, at 4:26 AM, John Reid <> wrote:
>
>> On 07/10/14 22:17, Bob Carpenter wrote:

Well I have the model as knitr output on a web-page: http://zemlinsky.mrc-bsu.cam.ac.uk/johns/Julia/GP.html so feel free to have a look. The code to generate data from the model is there together with how I'm using Stan to perform inference. The posterior plots show that the mean estimates of the pseudotime are not great but that the sample with the best lp__ does quite well. I'm happy to release the model under the MIT license. I hope you find it useful.

The MIT License (MIT)

Copyright (c) [2014] [John Reid]

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.