Request: Skew Student distribution

[Slawek Smyl]

Hi,

I am creating time-series models on some volatile (and some not) curves and found Student distribution really useful, as it can accommodate to the curve volatility.  I would like to experiment with simpler models, where trend or drift is expressed just with skewedness.

Currently I am doing it by collecting error distribution while fitting the model, and then in R utilizing this skew information for simulating future, but it is a inexact and a bit clumsy.

So, please add the Skew Student distribution - it would be very general and useful.

Regards,

Slawek

P.S. I read FAQ, and will try to do "increase lp__ trick". Still, my requests stands :-)

[Bob Carpenter]

Thanks for the note. I added a feature request to our issue tracker:

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

We're working on a lot of core components now and won't have
time to get to this any time soon. Having said that, it shouldn't
be hard to write given that we already have both student-t and skew
normal --- the one hitch may be that skew student-t would be more
parameters than any of our other distros and the testing framework
may not generalize.

As you found, the best bet is to write a function and use that:

functions {

real skew_student_t_log(real y, ...) {
real lp;
lp <- 0;
... increment lp ...
return lp;
}
}

Then you'll be able to use it in the language as

y ~ skew_student_t(...);

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

[Bob Carpenter]

Oh, and if you do define it, please add the definition to the issue
tracker (or mail back to the list if you don't have a GitHub account).
That way, we'll have a pattern to work from and other people will have
an easy to find plug-in function.

- Bob

> On Jun 8, 2015, at 3:37 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> Thanks for the note. I added a feature request to our issue tracker:
>
> https://github.com/stan-dev/stan/issues/1485
>
> We're working on a lot of core components now and won't have
> time to get to this any time soon. Having said that, it shouldn't
> be hard to write given that we already have both student-t and skew
> normal --- the one hitch may be that skew student-t would be more
> parameters than any of our other distros and the testing framework
> may not generalize.
>
> As you found, the best bet is to write a function and use that:
>
> functions {
>
> real skew_student_t_log(real y, ...) {
> real lp;
> lp <- 0;
> ... increment lp ...
> return lp;
> }
> }
>
> Then you'll be able to use it in the language as
>
> y ~ skew_student_t(...);
>
> - Bob
>
>
>> On Jun 8, 2015, at 11:37 AM, Slawek Smyl <slawe...@gmail.com> wrote:
>>

[Slawek Smyl]

Bob,

thank you very much for creating an issue.

I will (attempt to :-) ) create the skew_student_t function.  

Many thanks,

Slawek

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

[Slawek Smyl]

Hi Bob,

The function I implemented is as follows:

real skew_student_t_log(real y, real nu, real mu, real sigma, real skew) {

  real lp;

  lp <- 0;

  if (skew<=0)

    reject("skew has to be positive. Found skew=",skew);

  if (y<mu)

     lp<- student_t_log(y*skew, nu, mu*skew, sigma);

  else

     lp<- student_t_log(y/skew, nu, mu/skew, sigma);

  lp <- lp + log(2*skew)-log(1+skew*skew);

  return lp;

}

and I call it like this:

 

y ~ skew_student_t(nu, mu, sigma, skew);

 

It follows Fernandez, C. and Steel definition, also used by R package skewt (except that they do not use sigma or mu).

I did not attempt to code as per newer definition by Azzalini and Capitanio, because it requires multiplying the density by the cumulative distribution, and the latter has reportedly some numerical problems if used in this way: https://groups.google.com/forum/?fromgroups#!searchin/stan-users/skew$20student/stan-users/j44ipDdQep8/kIx4bC2GsA4J

 

That’s great, but e.g. I do not know how to efficiently sample from this distribution, and I would like to do it in generated quantities block.

So, looking forward for a full implementation of the Skew Student in Stan, yours :-)

Slawek

[Bob Carpenter]

Thanks. Your definition makes it clear how the skew's happening.

I love that you're testing inputs! You should also test that sigma >= 0.
Of course, if you're estimating skew, then you want to constrain it to
positive, just like sigma.

You could simplify the code a bit by assigning to lp straight off
and then directly returning in the conditionals. I even wrote
the tests to validate that all return paths end in return statements,
so you can't forget to return.

real skew_student_t_log(real y, real nu, real mu, real sigma, real skew) {
real lp;

if (skew <= 0)
reject("Skew has to be positive. Found skew=", skew);
if (sigma <= 0)
reject("Scale has to be positive. Found sigma=", sigma);
lp <- log(skew) - log1p(square(skew));
if (y < mu)
return lp + student_t_log(y * skew, nu, mu * skew, sigma);
else
return lp + student_t_log(y / skew, nu, mu / skew, sigma);
}

Not assigning to 0 and then adding will be more efficient for autodiff.

You can drop the multiplier 2 in the adjustment, because

log(2 * skew) = log(2) = log(skew)

and Stan only needs densities up to an additive constant.

And you can use the log1p() and square() operations for a bit more
efficiency and in the case of log1p(), robustness for small values of
skew (probably not a problem here if there's a reasonable prior on
skew or if it's fixed).

If you're only using this with constant values of skew, you can do
this:

real skew_student_t_log(real y, real nu, real mu, real sigma, real skew) {

if (skew <= 0)
reject("Skew has to be positive. Found skew=", skew);
if (sigma <= 0)
reject("Scale has to be positive. Found sigma=", sigma);
if (y < mu)
return student_t_log(y * skew, nu, mu * skew, sigma);
else
return student_t_log(y / skew, nu, mu / skew, sigma);
}

because the whole expression is just a constant. So far, there's no
way to test that in our functions (unlike in C++, where we do it all
the time).

And if sigma is a constant, you can probably get even more efficient
by scaling that rather than y. But that's getting carried away in
terms of clarity.

- Bob

[Aki Vehtari]

On Tuesday, June 9, 2015 at 8:59:38 AM UTC+3, Bob Carpenter wrote:

You can drop the multiplier 2 in the adjustment, because

  log(2 * skew) = log(2) = log(skew)

and Stan only needs densities up to an additive constant.  

 But you may need to include constants if you use skew_student_t_log in generated quantities block to compute terms for model comparison...

Aki

[Bob Carpenter]

Right.

If we want to compute WAIC, don't you want the likelihood on
the constrained scale, or is it OK to use the unconstrained scale?
If it's OK to use unconstrained (with Jacobians), then you should
be able to use lp__ for the joint unconstrained density if you don't
use any sampling statements (~), which automatically drop constants.

- Bob

[Slawek Smyl]

Hi,

I have just tried the newer definition of skew Student (by A Azzalini, R Arellano-Valle e.g. in  http://arxiv.org/pdf/1203.2376v1, eq.22) . It seems to be working, despite using the CDF, so maybe things improved since 2013 (see above in the thread) or I am using a nice range of inputs? The matching R package is called sn.

Also, I do not understand all the differences between the two formulas, but in my testing this one appears to be working better (the generated forecasts are better).

real skew_student_t_log(real y, real nu, real mu, real sigma, real skew) {

   real z; real zc;

   if (sigma <= 0)

     reject("Scale has to be positive. Found sigma=", sigma);

   z<- (y-mu)/sigma;

   zc<- skew*z*sqrt((nu+1)/(nu+square(z)));

   return -log(sigma) + student_t_log(z, nu, 0, 1) + student_t_cdf_log(zc, nu+1, 0, 1);

}

Regards,

Slawek

[Bob Carpenter]

> On Jun 10, 2015, at 12:42 PM, Slawek Smyl <slawe...@gmail.com> wrote:
>
> Hi,
> I have just tried the newer definition of skew Student (by A Azzalini, R Arellano-Valle e.g. in http://arxiv.org/pdf/1203.2376v1, eq.22) . It seems to be working, despite using the CDF, so maybe things improved since 2013 (see above in the thread) or I am using a nice range of inputs? The matching R package is called sn.
>
> Also, I do not understand all the differences between the two formulas,

The algebra's simpler --- it's the numerical properties of the
algorithms that'll get you. The key thing for efficiency is not to repeat computations
if you can avoid it. The main think for numerical accuracy is to
make sure that you use the compound functions wherever possible (e.g.,
log1p(x) instead of log(1 + x)), us the parameters with built-in transforms
wherever possible (e.g., bernoulli_logit(alpha) vs. bernoulli(inv_logit(alpha))
and always work on the log scale.

> but in my testing this one appears to be working better (the generated forecasts are better).

Unless there are numerical issues, you should get the same density
value and gradient, so how a function's defined will only matter if
it runs into numerical difficulties.

> real skew_student_t_log(real y, real nu, real mu, real sigma, real skew) {
>
> real z; real zc;
>
> if (sigma <= 0)
>
> reject("Scale has to be positive. Found sigma=", sigma);
>
>
> z<- (y-mu)/sigma;
>
> zc<- skew*z*sqrt((nu+1)/(nu+square(z)));
>
> return -log(sigma) + student_t_log(z, nu, 0, 1) + student_t_cdf_log(zc, nu+1, 0, 1);
>
> }

Thanks for sharing. That looks great. What happened to the test
for skew <= 0?

Now we can work on style :-) I'd add space between operators like *,
remove space between lines, and keep lines to 80 characters max. At
least those are the guidelines we try to follow ourselves in coding models
and in the C++ for Stan itself.

- Bob




- Bob

[Slawek Smyl]

Hi Bob,

In http://www.iew.uzh.ch/static/seminars/downloads/GHSkewTAasHaff05.pdf, in paragraph 3.2 the authors state
" There are several definitions presented in the literature that can be regarded as competing skew Student’s t-distributions" and they follow on showing that e.g. tail behavior is a bit different.
So they look to me like different functions, not just different expressions.

The skew<=0 was removed, because in the Azzalini's formulation, skew is symmetric around 0, and 0 gives you the normal Student distribution. In Fernandez and Steel's formula 1 is central.
I am not sure what is best translation between these two skews. 

Regards,
Slawek

- Bob

[Bob Carpenter]

Thanks --- I thought some of those looked like different functions!

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