Sign function?

[Stephen Martin]

Hey all,

Am I crazy, blind, or both? I couldn't find a sign() function in Stan. I rolled my own using x -> x/|x| , but I assumed Stan would have something similar. Is there not?

[Michael Betancourt]

Intentional — such functions break differentiability of the target

and the accuracy of any gradient-based algorithm.  You can

hack it, but we actively try to not make it easy given the really

bad consequences.

On Mar 31, 2017, at 6:45 PM, Stephen Martin <hwki...@gmail.com> wrote:

Hey all,

Am I crazy, blind, or both? I couldn't find a sign() function in Stan. I rolled my own using x -> x/|x| , but I assumed Stan would have something similar. Is there not?

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

[Andrew Gelman]

sign_x <- 2*step(x) - 1;

This should be cleaner than x/|x|.

[Stephen Martin]

Nevertheless, it's a useful function. Sure, it can break some gradients, but it doesn't break the gradient of the one I'm working on.
It's purely a computational oddity.  The problem is parameterizing a skew-normal in terms of mean, sigma, and gamma, rather than the direct parameterization of xi, omega, and alpha. I've worked out how to do it, and it performs very well, but one part of it needs a sign function.

parameters {
    real<lower = -1*gamma_max,upper=gamma_max> gamma;
    real<lower=0> sigma;
    vector[K] beta;
}
transformed parameters {
}
model {
    real r;
    real delta;
    real alpha;
    real omega;
    vector[N] xihat;

    r = gamma/fabs(gamma)*(2*fabs(gamma)/(4-pi()))^(1.0/3.0); <- Here. The non-sign version has (2*gamma/(4-pi))^(1/3), but this produces only positive gammas, because negative ones result in NaN. This is conceptually nonsense, so you just use the abs value, then put the sign back on. This is how the author of the SN papers reparameterized the model, to avoid a likelihood oddity at alpha=0, and to reduce correlation of parameters. In stan, this method works excellently, but there is no sign function. This is a useful case of the sign function, and it doesn't break the gradient; if anything, it actually fixes the gradient (without this method, you get NaN and 1000s of divergences. With the method, 0 errors, 0 divergences).
    delta = r/(b*sqrt(1 + r^2));
    alpha = delta/sqrt(1 - delta^2);
    omega = sigma/sqrt(1-b^2*delta^2);
    xihat = X*beta - omega*b*delta;

[Daniel Lee]

Hi Stephen, 

Having a sign function doesn't make this more convenient. It's just as easy to branch on gamma > 0. 

For this particular problem, have you tried the cbrt() function? I'm not sure if that behaves properly with negative numbers. 

Daniel

[Stephen Martin]

How does branching make this problem easier? It's more code than just throwing a fabs() and a sign function into the thing. And if gamma < 0, I'd still need to put its sign back into place, so the point is moot, no?

I haven't tried the cbrt function, but I'm not sure it'd make a difference.

Also, this is the exact parameterization used by Azzalini, it's well tested, and it performs well.
Attached are my derivations of the CP -> DP parameterization; heck, if the stan team wanted to encourage the use of the SN distribution, the CP parameterization is generally better and easier to understand. It's also got some benefits over the DP (stan's current method) when estimating models with some heterogeneity in skewness or variance; if you do regression on the location parameter (DP) with heterogeneity in variance or skewness, you /will/ get the 'wrong' answer, it'll be nonsensical. On CP this isn't so much an issue, because you're estimating means, not location parameters, and the model will be correct.  That is, conceptually, if you have a two-group scenario where the means are exactly equal, but they differ in variance, then the model will suggest the location parameters are very different (which is true, but probably irrelevant to your question, since location != expected value).

[Andrew Gelman]

I don't understand this discussion.  What's wrong with using 2*step(x) - 1?

<deriv.pdf>

[Stephen Martin]

Nothing. I was more taken back by the notion that the sign function is intentionally not implemented, even though there are instances where the sign function would be needed. It just seems odd to me; there are all sorts of ways people could break the continuity of the gradient using the functions already implemented, why disallow a legitimate mathematical function purely because it /could/ break continuity?

I actually like your solution better because I x/|x| isn't defined at 0, but your method is, I think.

[Bob Carpenter]

First, rather than taking signum(x) and then using it as
a multiplicand, it's much more effcient to branch in Stan
because branching doesn't add anything to the expression graph.
That's the real reason we don't have it. Branching on
parameters also breaks differentiability and we should be
raising warnings when users try it.

The issue we have with the built-ins is that we can put
warnings in the manual, but users won't see those. We could
have every use of such a function raise a warning and we
might in the future (like warnings on branching on parameters).

We have a bunch of step-like functions which break
differentiability and just have warnings next to them in
the manual. The functions we do supply are the ones that are
in the C++ standard library, and they have C++ signatures like:

double floor(double);

C++ allows double values to be assigned to integers (up to
user to deal with overflow), but Stan doesn't. That's the
built-in limitation in Stan and I did that intentionally so
as to preserve distinct typing.

The function we're seriously missing is the round to integer
function. Most examples we've seen of users who want to use
it would be breaking our samplers, so we're reluctant. We might
just build in some interpolation functions.

- Bob