limits, NaN, infinite parameters & variates

[Bob Carpenter]

Daniel and I were discussing what we want to do with
density and mass functions with limits. (I'm setting
aside CDFs until we get the prob functions sorted.)

We're still thrashing without a decision.

My own feeling is that we should throw errors early
and often and that every case below other than those
marked OK should just throw a domain error.

There seemed to be substantial sentiment the other way,
but I still don't really understand why. Maybe I confused
this with the CDF discussion, but I'd like to clarify now.
Maybe someone could be patient and re-explain to me why we
want to return -inf instead of raising errors.

Note that we've already given up trying to write a general
C++ probabability API as soon as we tied the density functions
to our auto-dif. So I think the only issue is how this plays
into Stan.

Here's the set of possible behaviors for a few probability
functions.


A. normal_log(y,mu,sigma);

sigma NaN error
sigma -inf error
sigma +inf error [1]
sigma (-inf,0] error [2]
sigma (0,+inf) OK

(y-mu) NaN error
(y-mu) -inf -inf
(y-mu) +inf -inf
(y-mu) fin OK

[1] Could limit sigma=inf case to be improper uniform
and return an arbitrary constant

[2] Could take sigma=0 to return +inf if y = mu
and -inf otherwise



B. beta_log(y,alpha,beta)

alpha NaN error
alpha [-inf,0) error
alpha 0 error
alpha (0,inf) OK
alpha inf -inf [4]

beta (just like alpha)

y (0,1) OK
y 0 -inf [4]
y 1 -inf [4]
y [-inf,0) error
y (1,inf) error

[4] alpha = inf, beta = finite, could lead to
+inf probability for y=1 and -inf otherwise;
similarly alpha=finite, beta=inf produces
the reverse


C. bernoulli_logit(y,alpha)

alpha NaN error
alpha -inf error [5]
alpha +inf error [5]
alpha finite OK

y < 0 error [6]
y > 1 error [6]
y {0,1} OK

[5] If alpha = -inf, could return +inf for y = 1 and -inf for y = 0

[6] Could return -inf for y not in {0,1} on the principle that
probability functions should be total for fixed valid params



Help!

- Bob

[Marcus Brubaker]

My general belief is that if the behaviour of a function can be consistently well defined then there's no need to throw an exception.  In most cases in stan there are no practical differences between the two cases for the domain of PDFs, but throwing exceptions reduces flexibility.  For instance, consider trying implement a mixture model where the components of the mixture have different supports.  E.G.,

p(x) = w_0 N(x|mu,sigma) + w_1 Beta(x|a,b) + w_3 Uniform(x|-1,0)

I'm not saying this is a common use case, but if a call to beta_log(x,a,b) throws an exception when x is outside of [0,1] then it would be hard (if not impossible) to implement it without a lot of extra work.  In contrast if beta_log returns -inf then it could be done in a relatively straightforward manner.

For distribution parameters I think exceptions make more sense because you could inadvertently get an improper prior in some cases and in many cases it's hard if not impossible to define things consistently.

Cheers,
Marcus

- Bob


--

[Matt Hoffman]

The reason we don't want to always throw an error is because infs can
happen naturally as a result of overflow during leapfrog steps. This
does indicate that the step size is too big, but that's something we
can and do deal with without dying.

I'm basically ok with Marcus's reasoning about not throwing domain
errors for random variables.

My feeling is that anywhere that a limit exists and makes sense we
should use it, since the typical case is that the simulation really
did "approach" 0 or infinity, and then got rounded off over the cliff.

The exception is where we'd get +inf density. We absolutely cannot
allow ourselves to get into such a state, because we'll never get out
again. I'm leaning towards saying we should actually return -inf when
we get a +inf. Which seems crazy, but it can be interpreted as
modifying the sampler to never transition to such states, which
actually preserves detailed balance (if not ergodicity).

A more principled reason is that in many cases you can justify it via
a change of variables to an unconstrained space. For example, consider
the gamma:
p(y) = exp{(a-1)*log(y) - b*y + a*log(b) - lgamma(a)}.
If you transform it so that x = log(y), then you get
p(x) = exp{a*x - b*exp(x) + a*log(b) - lgamma(a)}.
No more +infinities, and you get this limit behavior:
x=±inf (y=0 or inf): -inf.
a=0: -inf.
b=0: -inf.

My thoughts on specific distributions inline:

> A. normal_log(y,mu,sigma);

sigma < 0: error.
sigma = 0, y-mu = 0: error or -inf.
sigma = 0, y-mu ≠ 0: -inf.
sigma = +inf: -inf.
y-mu = ±inf: -inf.

The 2nd case gives infinite density, which is bad.
The 3rd case can arise organically as a result of underflow, so we
should allow the sampler to recover gracefully from it.
The 4th and 5th cases can arise as a result of overflow, so we should
allow the sampler to recover gracefully from them.

> B. beta_log(y,alpha,beta)

This is a tricky one. I say we go with the behavior on the
logit-transformed scale, i.e., think of x\in(-inf, inf) so that
p(x) = exp(alpha*log(x) - beta*log(1-x) - lgamma(alpha) - lgamma(beta)
+ lgamma(alpha + beta)).
Then things are clearer:
x=±inf (y=0 or 1): -inf
alpha or beta = 0: -inf
Errors if alpha or beta < 0.

> C. bernoulli_logit(y,alpha)

This is a PMF, not a density, so numerical issues are less of a big
deal. Edge cases:
alpha = inf, y=1: 0
alpha = -inf, y=1: -inf
alpha = -inf, y=0: 0
alpha = inf, y=0: -inf

Matt

[Ben Goodrich]

On Thursday, October 18, 2012 9:28:58 PM UTC-4, Matt Hoffman wrote:

> B.  beta_log(y,alpha,beta)

This is a tricky one. I say we go with the behavior on the
logit-transformed scale, i.e., think of x\in(-inf, inf) so that
p(x) = exp(alpha*log(x) - beta*log(1-x) - lgamma(alpha) - lgamma(beta)
+ lgamma(alpha + beta)).
Then things are clearer:
x=±inf (y=0 or 1): -inf
alpha or beta = 0: -inf
Errors if alpha or beta < 0.

 

As both alpha and beta approach 0, you can define the beta as a Bernoulli with p=0.5. Otherwise, I would agree.

Ben

[Matt Hoffman]

Right, the only thing I worry about in that case is if they aren't
approaching 0 at the same rate (which we have no way of tracking), in
which case you wind up with a different Bernoulli. Anyway, we'd still
have to worry about how to cope with infinite density.

> --
>
>

[Daniel Lee]

On Oct 18, 2012, at 9:28 PM, Matt Hoffman <mdho...@cs.princeton.edu> wrote:
>> B. beta_log(y,alpha,beta)
>
> This is a tricky one. I say we go with the behavior on the
> logit-transformed scale, i.e., think of x\in(-inf, inf) so that
> p(x) = exp(alpha*log(x) - beta*log(1-x) - lgamma(alpha) - lgamma(beta)
> + lgamma(alpha + beta)).
> Then things are clearer:
> x=±inf (y=0 or 1): -inf
> alpha or beta = 0: -inf
> Errors if alpha or beta < 0.

Matt, you skirted a point of difficulty that Bob and I had when thinking about this. What should the behavior be when y < 0 or y > 1? Both -inf and error make sense. I lean towards error because values of y outside (0, 1) indicates that something was misspecified.

[Matt Hoffman]

That's been my thought too, although Marcus brings up at least one
case (mixture modeling) where -inf would make sense. But it's possible
that the utility of that sort of case is small compared to the utility
of throwing an interpretable error message in the more common case
where someone has misspecified something.

Matt

> --
>
>

[Marcus Brubaker]

Once we get the code in place to give model file line numbers with exceptions, we will also be able to easily check for lp__ == -inf and issue meaningful error messages without needing the *_log functions to throw exceptions.  We will still need to make that check due to potential underflow or programming issues with people directly manipulating lp__.

The other possibility is to make the behaviour configurable, eg, with a template parameter.  Calls that happen as the result of a ~ statement could throw exceptions and all other calls return -inf for out of bounds behaviour.

Cheers,
Marcus

--

[Matt Hoffman]

On Fri, Oct 19, 2012 at 8:31 AM, Marcus Brubaker
<marcus.a...@gmail.com> wrote:
> Once we get the code in place to give model file line numbers with
> exceptions, we will also be able to easily check for lp__ == -inf and issue
> meaningful error messages without needing the *_log functions to throw
> exceptions. We will still need to make that check due to potential
> underflow or programming issues with people directly manipulating lp__.

I don't think we shouldn't be throwing exceptions (line numbers or no)
just because lp__ == -inf, for the reasons I mentioned above. In my
opinion, we should only be throwing exceptions when we know the
problem is due to a modeling problem or due to something that could
only be due to a bug in Stan, not due to numerical issues. We really
need to be able to gracefully handle the scenario where a parameter
gets sufficiently large or small that lp__ overflows to -inf. Or, for
that matter, we'd like to potentially be able to support posteriors
with holes in them.

> The other possibility is to make the behaviour configurable, eg, with a
> template parameter. Calls that happen as the result of a ~ statement could
> throw exceptions and all other calls return -inf for out of bounds
> behaviour.

I like this idea.

Matt

[Bob Carpenter]

On 10/19/12 11:38 AM, Matt Hoffman wrote:
...

> I don't think we shouldn't be throwing exceptions (line numbers or no)
> just because lp__ == -inf, for the reasons I mentioned above. In my
> opinion, we should only be throwing exceptions when we know the
> problem is due to a modeling problem or due to something that could
> only be due to a bug in Stan, not due to numerical issues. We really
> need to be able to gracefully handle the scenario where a parameter
> gets sufficiently large or small that lp__ overflows to -inf.

By graceful here, do you just mean a normal (non-exception)
return from the function?

I keep thinking of the behavior in HMC, which is essentially
the same either way because we Metropolis reject -inf states.

We can already configure that at the Boost error level by having
domain errors not throw exceptions.

> Or, for
> that matter, we'd like to potentially be able to support posteriors
> with holes in them.

You mean in the above sense of just having a non-exceptional
return?

>> The other possibility is to make the behaviour configurable, eg, with a
>> template parameter. Calls that happen as the result of a ~ statement could
>> throw exceptions and all other calls return -inf for out of bounds
>> behaviour.
>
> I like this idea.

I think what Daniel may be suggesting in the way of more
configurability is to have a separate kind of Boost error
behavior, say "neg_inf_result", which we raise when there is
a -inf result. The policy could then configure whether to
proceed or raise an exception before returning -inf.

All else being equal, I like configurability, too. The main
issue is the growing complexity of the templates, particularly
how we have to hack pre-C++11 function template defaults.

(Speaking of which, one way to fix those would be to move
to static class functions rather than standalone functions.)

- Bob

[Matt Hoffman]

Sorry, I've been conflating exceptions and fatal errors. I'm ok with
non-fatal exceptions, and defer to others' opinions on whether it's
better to sound an alert when we hit -inf or not. My strong feelings
were just about fatal exceptions—sorry for muddying the discussion.

Matt

> --
>
>

[Bob Carpenter]

n 10/19/12 2:24 PM, Matt Hoffman wrote:
> Sorry, I've been conflating exceptions and fatal errors. I'm ok with
> non-fatal exceptions, and defer to others' opinions on whether it's
> better to sound an alert when we hit -inf or not. My strong feelings

> were just about fatal exceptions�sorry for muddying the discussion.

Hang on! You made me (and Daniel) think more deeply
about all of this.

I'd be happy going totally the other way with errors
and instead follow the std libs and

(a) return NaN when we get a NaN input

(b) return the appropriate finite, -inf, +inf or NaN depending
on the input configs

(c) generalizing functions as much as possible [e.g., symmetrizing
infinities on y and mu in normal(y|mu,sigma)].

As Daniel said, this is WAY simpler than what we have in
place now. It'd cut out whole files full of error handling code
and "redundant" function signatures for default error handling.

We'd lose those "non-fatal warning (error: found NaN in ...)
and so on. But those were just confusing people and making
them apprehensive about whether Stan was working.

Sampling output should be exactly the same (HMC/NUTS will reject
an -inf or NaN state). I don't think there's an efficiency
issue because this is all in "error" handling.

- Bob

[Bob Carpenter]

On 10/18/12 9:08 PM, Marcus Brubaker wrote:
> My general belief is that if the behaviour of a function can be consistently well defined then there's no need to throw
> an exception. In most cases in stan there are no practical differences between the two cases for the domain of PDFs,
> but throwing exceptions reduces flexibility.

Let's say it's differently abled flexibility-wise.

The exception makes it easier to report an edge case
likely arising due to an error. In general, exceptions
are nicely flexible if you want to throw from multiple
functions and have a single catcher higher up.

In effect, our rejection for -inf or NaN values has the
same effect as the catch clause higher up, so the
exceptions aren't buying us much other than error reporting
and early termination (though I haven't looked into how.

The boundary return makes it easier to do arithmetic.
I think the case you give next is a good example!

Returning NaN or -inf should be at least as fast as
try/catch -- I started reading up on try/catch overhead
and quickly got in deeper than I wanted to. I think
the answer is that if the other ops are expensive, try/catch
isn't that expensive if there's no catch, but it is some overhead
in tracking state.

> For instance, consider trying implement a mixture model where the
> components of the mixture have different supports. E.G.,
>
> p(x) = w_0 N(x|mu,sigma) + w_1 Beta(x|a,b) + w_3 Uniform(x|-1,0)
>
> I'm not saying this is a common use case, but if a call to beta_log(x,a,b) throws an exception when x is outside of
> [0,1] then it would be hard (if not impossible) to implement it without a lot of extra work. In contrast if beta_log
> returns -inf then it could be done in a relatively straightforward manner.

This is a good example of where you can go to infinity and
back again:

log_sum_exp( log(w_0) + normal_log(x,mu,sigma),
log(w_1) + beta_log(x,a,b),
log(w_3) + uniform_log(x|-1,0) );

If any of the terms to log_sum_exp is -inf, it just gets
dropped, becase exp(-inf) = 0.

> For distribution parameters I think exceptions make more sense because you could inadvertently get an improper prior in
> some cases and in many cases it's hard if not impossible to define things consistently.

This is where Matt was arguing before we should have symmetry
between the errors thrown by:

normal_log(y, mu, sigma)

normal_log(0, y - mu, sigma);

normal_log(0, (y - mu)/sigma, 1);

- Bob