Issue 76 in stan: beta(0.8,0.8) nearly always generates ~0

[st...@googlecode.com]

Status: New
Owner: ----
Labels: Type-Defect Priority-Medium

New issue 76 by lewf...@gmail.com: beta(0.8,0.8) nearly always generates ~0
http://code.google.com/p/stan/issues/detail?id=76

What steps will reproduce the problem?

1) cd <stan-home>
2) make ../stan-test/beta-0.8
3) cd ../stan-test
4) ./beta-0.8 --samples="samples-0.8.csv"

What is the expected output? What do you see instead?

I expect the values for theta in samples-0.8.csv to be distributed
according to beta(0.8,0.8), approximately uniformly between 0 and 1, but
instead they are all tiny values close to zero (around 1e-16). This is
unexpected because this should be a very easy model to sample from.

Some other observations:

-Initializing theta to be 0.99 doesn't help.
-On a small fraction of runs, theta will be a reasonable value for a few
steps and then goes to a tiny value close to zero and stays there.
-When changing the model slightly to beta(1.0,1.0) the problem does not
seem to occur.

What version of the product are you using? On what operating system?
Stan 1.0.2
g++ (Ubuntu/Linaro 4.6.1-9ubuntu3) 4.6.1
Ubuntu 11.10
i686 32-bits
2gb RAM

Attachments:
beta-0.8.stan 88 bytes
samples-0.8.csv 29.6 KB

[st...@googlecode.com]

Updates:
Status: Accepted
Owner: bea...@alum.mit.edu
Cc: ca...@alias-i.com

Comment #1 on issue 76 by bea...@alum.mit.edu: beta(0.8,0.8) nearly always
generates ~0
http://code.google.com/p/stan/issues/detail?id=76

seed=2605059476 is problematic for this model. We'll have to figure out
what's going on.

There isn't anything wrong with the model. For example, seed=2614914497
(randomly generated) works fine.

If we aren't sampling properly, we should at least be throwing up messages.

[st...@googlecode.com]

Comment #2 on issue 76 by goodrich...@gmail.com: beta(0.8,0.8) nearly

always generates ~0
http://code.google.com/p/stan/issues/detail?id=76

It appears to be poor numeric behavior at the endpoints of the [0,1]
interval. For example, this model

transformed data {
real<lower=0> eps;
eps <- 0.0000001;
}
parameters {
real<lower=eps,upper=1-eps> theta;
}
model {
theta ~ beta(0.8,0.8);
}

seems to be fine with --seed=2605059476

goodrich@CYBERPOWERPC:/tmp/SU/beta$ R --quiet
-e "summary(read.csv('samples.csv', comment.char='#'))"
> summary(read.csv('samples.csv', comment.char='#'))
lp__ treedepth__ stepsize__ theta
Min. :-6.280 Min. :1.000 Min. :4.463 Min. :0.003108
1st Qu.:-2.070 1st Qu.:1.000 1st Qu.:4.463 1st Qu.:0.186391
Median :-1.488 Median :1.000 Median :4.463 Median :0.495658
Mean :-1.774 Mean :1.179 Mean :4.463 Mean :0.501181
3rd Qu.:-1.196 3rd Qu.:1.000 3rd Qu.:4.463 3rd Qu.:0.802087
Max. :-1.109 Max. :2.000 Max. :4.463 Max. :0.999610

[Ben Goodrich]

Which raises the general question, if a user writes

real<lower=a,upper=b> theta;

should that silently parse in C++  to a + guard and / or b - guard for some extremely small (perhaps tweakable) value of guard?

Ben

[Daniel Lee]

Although that may solve the issue, I think something weird is going
on. Why are we initializing at 0 if the inits are on a -2 to +2 on the
transformed scale? Shouldn't that put the initial point somewhere
between invlogit(-2) and invlogit(2)? or 0.12 and 0.88?

> --
>
>

[st...@googlecode.com]

Updates:
Status: Fixed
Cc: -ca...@alias-i.com

Comment #3 on issue 76 by bea...@alum.mit.edu: beta(0.8,0.8) nearly always
generates ~0
http://code.google.com/p/stan/issues/detail?id=76

We updated the transform code to bound away from zero. This should solve
the problem for this model and others using the lub_transform() code in
prob/transform.hpp

[Ben Goodrich]

You can go with a lot smaller of a guard at zero than e-15. In R,

dbeta(1e-323,.8,.8,log=TRUE) # 148.3327

but anything smaller overflows to Inf.

Ben

[Daniel Lee]

It used to be 1e-100 before I changed it.

> --
>
>

[Ben Goodrich]

We must be losing the precision somewhere in the transformation then. Anyway, we can worry about it later.

[Daniel Lee]

It might not be the transformation, it could be the derivative.

> --
>
>

[Bob Carpenter]

On 10/25/12 3:18 PM, Ben Goodrich wrote:> You can go with a lot smaller of a guard at zero than e-15. In R,

>
> dbeta(1e-323,.8,.8,log=TRUE) # 148.3327
>
> but anything smaller overflows to Inf.

It didn't work at the original 1e-100 bound.

When Matt originally (I think it was Matt) put those
bounds in, he used 1 - 1e-15 for the upper bound and
1e-100 for the lower bound. The upper bound has to
be tighter because you can get closer to 0 with floating
point numbers than to 1.

The beta(x|0.8,0.8) turns into a function f(x) s.t.

f(x) = x^(0.8 - 1) * (1 - x)^(0.8 - 1)

log f(x) = -0.2 * (log x + log (1 - x))

d/dx log f(x) = -0.2 * ( 1 / x - 1 / (1 - x))

If x is very tiny or very close 1, big problems (pun
intended).

- Bob

[Ben Goodrich]

I figured something along those lines. But R is evaluating the same function with more robustness, albeit without the transformation or the derivative. Also, at some point we should special case the beta(theta,1,1) call.

Ben

[Bob Carpenter]

Why do you say that? In Stan, we get:

beta_log(1e-100,0.8,0.8) = 45.635
beta_log(1e-300,0.8,0.8) = 137.738
beta_log(1e-323,0.8,0.8) = 148.333

It's the same to 6 decimal places, at least!

It's the gradient that's problematic, not the value.

> Also, at some point we should special case the beta(theta,1,1) call.

In this case, we'll evaluate

0 * log(x) + 0 * log(1 - x)

and add a useless 0 term to the auto-dif
stack for x, I think.

- Bob

[Ben Goodrich]

That's what I mean, we can't get the full precision of the function because we have to calculate the gradient.
 

 > Also, at some point we should special case the beta(theta,1,1) call.

In this case, we'll evaluate

   0 * log(x) + 0 * log(1 - x)

and add a useless 0 term to the auto-dif
stack for x, I think.

 
But if either log(x) or log(1 - x) blows up to -inf, then 0 * -inf will go nan, when the value is well-defined to be 0.0.

Ben

[Bob Carpenter]

On 10/25/12 4:26 PM, Ben Goodrich wrote:
> On Thursday, October 25, 2012 4:19:07 PM UTC-4, Bob Carpenter wrote:
>
> On 10/25/12 4:08 PM, Ben Goodrich wrote:
> > On Thursday, October 25, 2012 3:53:56 PM UTC-4, Bob Carpenter wrote:
> >
> > On 10/25/12 3:18 PM, Ben Goodrich wrote:> You can go with a lot smaller of a guard at zero than e-15. In R,

> > I figured something along those lines. But R is evaluating the same function with more robustness, albeit without
> the
> > transformation or the derivative.
>
> Why do you say that? In Stan, we get:
>
> beta_log(1e-100,0.8,0.8) = 45.635
> beta_log(1e-300,0.8,0.8) = 137.738
> beta_log(1e-323,0.8,0.8) = 148.333
>
> It's the same to 6 decimal places, at least!
>
> It's the gradient that's problematic, not the value.
>
>
> That's what I mean, we can't get the full precision of the function because we have to calculate the gradient.

Stan evaluates the log(beta(x|alpha,beta)) function just fine.
If you provide double-valued inputs, no gradient gets calculated,
and the result is the same as R's (it's a simple function, after
all).

If you use var inputs, you get the right value of the function,
but the gradients will overflow and HMC gets hosed.

Is there something you think we can make better here?

- Bob

[Marcus Brubaker]

In those cases you can replace the true gradient with one which approximates the true gradient but is numerically stable.  A simple case of this would be to put an upper limit on the gradient magnitude.

Cheers,

Marcus

[Ben Goodrich]

 
Something must be going on with the combination of the transformation and the gradient. Even for very small theta, the derivative with respect to theta should be decently behaved. Like

> theta <- 1e-300
> -0.2 * (1 / theta - 1 / (1 - theta))
[1] -2e+299

The 1 - x piece underflows to 1.0 but that is okay. I guess that when you multiply that by the derivative of theta wrt the unbounded parameter, there could be problems. But even that looks like it should work

> unbounded <- -700
> theta <- plogis(unbounded)
> -0.2 * (1 / theta - 1 / (1 - theta)) * dlogis(unbounded)
[1] -0.2

Ben