limit behavior of cdfs

[Bob Carpenter]

I just added tests for limit behavior
for cdfs:

OK: normal_cdf(-inf,0,1) = 0
OK: normal_cdf(+inf,0,1) = 1

OK: exponential_cdf(0,1.5) = 0
OK: exponential_cdf(+inf,1.5) = 1

OK: lognormal_cdf(0,0,1) = 0
OK: lognormal_cdf(+inf,0,1) = 1

OK: weibull_cdf(0,1,1) = 0
OK: weibull_cdf(+inf,1,1) = 1

(I had to fix +inf for weibull_cdf by removing
the finiteness test on the variate y.)

The beyond-the-limit-of-support items fail
now for two of the three cdfs:

FAIL: exponential_cdf(-inf,1.5) = 0
FAIL: lognormal_cdf(-inf,0,1) = 0
OK: weibull_cdf(-inf,1,1) = 0

I don't care which way the decision goes, but
I do think they should all behave the same way.

- Bob

[Marcus Brubaker]

Those should all pass for sure as an exponential distribution can be
considered to be defined on the entire real line, but with zero
support for values less than zero. In fact, they should be == 0 for
any point less than 0.

Cheers,
Marcus

[Bob Carpenter]

Andrew -- do you care about this one way or the other?

We need to decide once and for all what happens for
variates outside their support for cdfs, pdfs, and pmfs.

The current behavior for probability functions with
variates outside their support is to raise a domain error.
For example, calling beta_log(-1,1,1), where -1 is
outside the support of [0,1], raises a std::domain_error.
The error effectively cancels any remaining computation
and rejects the step in the leapfrog algorithm.

R, on the other hand, just returns 0, as in:

> dbeta(0.5,2,2)
[1] 1.5
> dbeta(-1,2,2)
[1] 0
> dbeta(1,2,2)
[1] 0

I very strongly feel that we should do the same thing
for cdfs as for pdfs and pmfs.

Marcus came down on the side of returning 0 (-infinity on
the log scale). I'm happy to push this change through
all of our functions, but want more of a consensus before
undertaking the change.

I was operating under the catch-errors-early principle
where as soon as we find a coding error, we raise an
exception, halt computation, and alert the user with an
informative message.

So what I'd like to hear from those who think we should
just return 0 is when a user would intend to provide
a probability or distribution function with a variate outside
of its support.

I'm also willing to just go along with the consensus minus
me if everyone's too busy to respond but feels strongly
about it.

Previous discussion below.

- Bob

[Marcus Brubaker]

Just to expand on this point a little I view CDFs and PDFs as somewhat
different. A CDF can be mathematically well defined for any value on
the real line and while it might seem silly to write x ~
exponential(lambda) T[-100,] I can't see a reason to disallow it. For
instance, if you had a lower bound parameter b so

b ~ normal(0,1)
x ~ exponential(lambda) T[b,]

where b is used elsewhere in the model it seems unnecessary to cause
an error either forcing the user to construct a variable which is just
max(b,0) or slowing the sampler down due to a lot of rejected
trajectories.

With PDFs there are two cases I can think of. One is the typical use-case

model {
x ~ exponential(lambda)
}

where x is data or some other parameter. If x < 0 then
exponential_log will return either -infinity or throw an exception.
The exception will get caught and NUTS will terminate that end of the
trajectory the -infinity will propagate to lp__ and cause NUTS to
similarly terminate that end of the trajectory. The result of either
option is the same, but the exception will result in some small amount
of saved computation on the step because everything after that point
won't be evaluated.

The flip side of things is that someone may want to use the density
functions on their own. For instance, a model I played with recently
had an exponential decay after an event at time t0 which could be well
modeled with exponential_pdf(t - t0,lambda) which relies on the
behaviour being proper when the variate is less than zero.

I'm not really certain what the right answer is here for the PDFs as
both have some pros and cons. (I think it's clearer with the CDFs.)
A third option would be to have a parameter or different functions for
different behaviour but that requires more work.

Cheers,
Marcus

> --
>
>

[Andrew Gelman]

I like the idea of returning 0.

> --
>

[Ben Goodrich]

I agree with Marcus.

Ben

[Bob Carpenter]

I think the pdfs just being derivatives of
the cdfs is a strong argument for consistency.

For instance, if exp_cdf(x|theta) is the exponential
distribution's cdf, then if we think of it as being
defined on (-inf,inf), then the pdf should be defined
wherever the cdf is differentiable. And we've just
been filling in natural points when things aren't
differentiable.

- Bob

[Bob Carpenter]

OK -- sounds like Ben, Andrew and Marcus
all agree, which is enough of a consensus
for me (Matt's usually gone that way on these
issues, too). I'm leaning the other way because
I still haven't seen a motivating example.

n the fence, and I don't
think Daniel cared one way or the other.

Daniel and I will get busy changing all the
existing code and doc and the change will go
out in 1.0.2 or 1.0.3.

Is this worth dealing 1.0.2 a day or two for?

Now is more consistent, later lets us get
the patches in 1.0.2 out more quickly.

I'm leaning toward later, as I want to cut down
on the number of versions we push out at the
expense of leaving buggy ones out there.
I could be convinced to lean the other way
if we could automate/coordinate the releases
on the RStan side a little more easily.

- Bob

[Marcus Brubaker]

I'm in the release early, release often camp myself. Put out 1.0.2
ASAP since it resolves some bigger issues and then do this overhaul
for 1.0.3. It's not like we have to pay for version numbers or
anything ;)

Cheers,
Marcus

> --
>
>

[Bob Carpenter]

On 9/24/12 5:56 PM, Marcus Brubaker wrote:
> I'm in the release early, release often camp myself. Put out 1.0.2
> ASAP since it resolves some bigger issues and then do this overhaul
> for 1.0.3. It's not like we have to pay for version numbers or
> anything ;)

We don't have to pay in dollars, but we do pay
in user effort for all the installs and also
the fact that it makes the project look a bit
hysterical.

And we pay in my effort at making sure the doc's
OK, changing version numbers in the code, writing
the release notes, e-mailing the lists, actually doing the
upload. And that's just Stan. For RStan, there's also
Jiqiang's effort along the same lines.

- Bob