Saving samples of quantities computed based on variables defined in the model block

[Paul Buerkner]

Hi all,

I am currently implementing the generalized extreme value distribution in brms, which proved to be rather complicated.
Basically, I have to constrain the shape parameter xi to a certain interval based on other (predicted) model parameters to avoid
leaving the definition space of the response variable, which may cause many divergent transitions (the .stan file is attached).

These other predicted parameters are defined and computed in the model block, and not in transformed parameters,
to avoid that rstan is writing them to a csv-file in the temp folder even if not returned to the user later on
(at least this was rstan's behavior a few months ago when I moved most stuff to the model block).
These csv-files will get really big for huge models, which reportedly caused problems.

Long story short, I want to save the samples of the shape parameter xi, even though it depends on variables
defined in the model block. Or, I want to find a way to make sure that samples of quantities, which are never returned
to the user anyway, are not written to csv-files in the temp folder when defined in transformed parameters.

Many thanks for your help!

Best,
Paul

[Andrew Gelman]

Hi, my quick question (without having seen any of the models) is that you might want stronger priors.  Many times when a model is hard to fit, the problem is that the fitting algorithm is spending time in weird, remote areas of the posterior that correspond to data distributions which are not of any substantive interest.  Even a weak prior can be enough to make a big difference.

A

--
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.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<gev2.stan>

[Paul Buerkner]

I have tried informative priors, but that doesn't stop Stan from leaving the definition space. The GEV distribution is defined on y values satisfying

1 + xi (y - mu) / sigma > 0

which kind of sucks since it depends on all three parameters of the distribution.

[Jonah Gabry]

Hey Paul, one way to save xi would be to write a function that computes xi (and the intermediate stuff you need to eventually compute xi) and then use that function in the model block and then again in the generated quantities block, where you can create the xi that's actually saved. This will duplicate some computations but I'm not sure if there's a way around that at the moment (although maybe someone else will have a brilliant work-around). 

Jonah

[Bob Carpenter]

Definitely what Jonah said---recomputing in generated quantities
is very cheap because there are no derivatives and it's only
done once per iteration.

There's no way to get around that scaling by hand that you're
doing because we don't yet allow varying upper and lower bounds
on a variable (on our to-do list). You sidestepped the Jacobian
calc by putting the prior on the unconstrained form of xi. I
don't know if that makes sense in your generative model, but you
could put the prior directly on xi if you include the Jacobian
(which is described in the manual---we mean to expose those to
the language soon, too, for just this kind of application).

- Bob

[andre.p...@googlemail.com]

Paul,

I would fit towards the mean not towards the parameters of your distribution. Since the mean for GEV is not defined for zeta >= 1(shape)

maybe go for the mode instead. This doesn't change things for correctness point of view, but from sampling.

The inv_logit with max/min, doesn't feel good to me in MCMC. Did you try to constrain your parameters? So they get truncated.

Andre

[Bob Carpenter]

Presumably the bad feeling is from the discontinuities
in min and max? That does seem like it could be a problem
in the constraining transform:

real scale_xi(real xi, vector y, vector mu, vector sigma) {
vector[2] bounds;
real lb;
real ub;
bounds[1] = - inv(min((y - mu) ./ sigma));
bounds[2] = - inv(max((y - mu) ./ sigma));
lb = min(bounds);
ub = max(bounds);
return inv_logit(xi) * (ub - lb) + lb;
}

- Bob

[Paul Buerkner]

Thanks for all those comments. The problem I see is that when using the generated quantities block, we again get everything written to csv-files (even if not actually saved afterwards) or am I mistaken?

[Paul Buerkner]

I am not sure whether this is a problem. With min max I am just finding the range of values that are allowed for xi, so that the definition space of y is valid for all data thus getting rid of many divergent transitions..

[Álvaro M]

Hi.

If you define a variable  in a local block in generated quantities with curly bracket, that variable is not save.

[Jonah Gabry]

On Friday, February 3, 2017 at 3:34:54 AM UTC-5, Paul Buerkner wrote:

Thanks for all those comments. The problem I see is that when using the generated quantities block, we again get everything written to csv-files (even if not actually saved afterwards) or am I mistaken?

None of the intermediate quantities created inside the function you write would get written.

[Bob Carpenter]

Suppose you have a 2-vector v and the function y = min(v).
Then the derivatives dy/dv[1] and dy/dv[2] are discontinuous
at the point where v[1] = v[2]. So you have to cross a discontinuity
in derivatives to move from (0.99, 1) to (1.01, 1). That can
be problematic. It's not that HMC will do the wrong thing, just
that it'll revert to a rejection sampler if it can't smoothly
cross boundaries; how bad that is will vary from problem to problem.

- Bob


> On Feb 3, 2017, at 3:39 AM, Paul Buerkner <paul.b...@gmail.com> wrote:
>
> I am not sure that this is a problem. With min max I am just finding the range of values that are allowed for xi, so that the definition space of y is valid for all data thus getting rid of many divergent transitions..

[Paul Buerkner]

Thanks again, these comments were super helpful!

[Bob Carpenter]

Glad it's useful.

I learn a lot from our lists, so I try to contribute
when I know the answer. I used to always try to answer,
but have scaled way back now that many more people are
jumping in.

- Bob

[Paul Buerkner]

I guess, I will have stick to this implementation for now, unless someone comes up with a better solution avoiding these discontinuities.

- Paul