Re: [stan-dev] numerical stability tips

[Bob Carpenter]

[move to stan-users; bcc-ing Sharad]

In general, the best thing to do is work on the log scale and
always use CCDFs where available. Unfortunately, in Stan,
our log CDFs and log CCDFs are implemented naively, so that's
not going to help.

The problem is that you can represent numbers very close to 0
using IEEE floating point, but not numbers very close to 1. That's
why we get both erf() and erfc() and log() and log1p().

Aside from underflow/overflow, which the erfc() and log1p()s of the world
try to mitigate, another big problem for arithmetic is subtraction of numbers
that are close to each other or on different scales. For instance, if beta_cdf
is near 1, then (1 - beta_cdf(...)) is going to cause what's known
as "catastrophic cancellation" --- you lose precision equal to the
shared leading digits.

The only course of action for a complex expression is to rearrange
terms so it's more stable and work on the log scale for as long as
possible, only converting back at the very end. You might find
that everything cancels nicely given the a and a+1 in the beta.

You can test overflow and underflow over input ranges of interest.
Accuracy requires a higher-precision package (e.g., Boost uses
something like 1000 decimal place accuracy arithmetic in some of
their tess).

You might find this form more stable:

exp(log(1 - beta_cdf(x, a+1, b))
- log(1 - beta_cdf(x, a, b))
+ log(a)
- log(a + b))

=

exp(beta_ccdf_log(x, a+1, b)
beta_ccdf_log(x, a, b)
+ log(a)
- log(a + b))

What we could really use is someone to implement high-precision
log CDFs and CCDFs.

- Bob

> On Aug 11, 2015, at 11:14 AM, Alp Kucukelbir <akucu...@gmail.com> wrote:
>
> hi folks,
>
> i've been communicating with sharad goel, who wants to use the function below in a Stan program.
>
> functions {
> // mean of beta(a,b) conditioned on being larger than x
> real beta_conditional_mean(real x, real a, real b) {
> return (1 - beta_cdf(x, a+1, b)) / (1 - beta_cdf(x, a, b)) * a / (a+b);
> }
> }
>
> any tips to make this more numerically stable? i think i'm mainly worried about when a and b are small, but i may be missing something here.
>
> cheers
> alp
>
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Michael Betancourt]

In one of the recent pull requests to the math library
I updated beta_cdf to use a _far_ more stable implementation,
so make sure to update to the most reason version.
But you’ll still have the cancellation problems that
Bob noted.

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

[Luc Coffeng]

+1! (I won't be much help though)

On Tuesday, 11 August 2015 19:05:57 UTC+2, Bob Carpenter wrote:

[move to stan-users; bcc-ing Sharad]