Fwd: [stan-users] beta distribution

[Ben Goodrich]

We really have to do something about the numerical stability of the beta distribution. It even blows up when you restrict the shape parameters to (0.01, 10).

Ben

On Sunday, February 24, 2013 3:04:50 AM UTC-5, Ben Goodrich wrote:

On Sunday, February 24, 2013 2:57:47 AM UTC-5, Bob Carpenter wrote:

I don't understand why you're making the priors for 
a and b so complicated. 

There is that. But I can't get any beta model to work on these data.

goodrich@CYBERPOWERPC:/tmp/SU/rubi$ ./test --test_grad --data=dat.txt --init=0

TEST GRADIENT MODE

 Log probability=-nan

 param idx           value           model     finite diff           error

         0               0            -inf            -inf            -nan

         1               0         15.8735            -nan            -nan

goodrich@CYBERPOWERPC:/tmp/SU/rubi$ ./test --test_grad --data=dat.txt

Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

 

See attachments.

[Matt Hoffman]

One problem is that we translate from unconstrained to constrained to
evaluate the log-density. This shows up with the gamma(0.01, 1) too,
whose left tail looks almost identical to the beta(0.01, 1)
distribution. If we have
x ~ gamma(0.01, 1)
then a typical draw for x is in the range of exp(-1000) to 5 or so.
But of course exp(-1000)=0 at double precision, and log(0) is outside
of the support of the unconstrained variable.

A workaround is to do something like this for the gamma(a, b):

parameters {
real logx;
}
transformed parameters {
real x;
x <- exp(logx);
}
model {
lp__ <- lp__ + a * logx - b * x;
}

or this for the beta(a, b):

parameters {
real logitx;
}
transformed parameters {
real x;
x <- inv_logit(logitx);
}
model {
lp__ <- lp__ - a * log(1 + exp(-logitx)) - b * log(1 + exp(logitx));
}

But that's not ideal. I'm not sure what we can do about this within
our current architecture, which really wants to transform beta r.v.s
to [0, 1] and gamma r.v.s to [0, inf].

Matt

> --
> 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/groups/opt_out.
>
>

[Ben Goodrich]

On Monday, February 25, 2013 4:57:42 PM UTC-5, Matt Hoffman wrote:

But that's not ideal. I'm not sure what we can do about this within
our current architecture, which really wants to transform beta r.v.s
to [0, 1] and gamma r.v.s to [0, inf].

All true, but in this particular case, the RV was data and not being transformed. So, it seems that for alpha very small and beta very large and x close to zero, the gradient of the log density blows up. Possibly the beta distribution is inappropriate for these data but it would be nice if we could come up with a solution to our numerical problems.

Ben

[Matt Hoffman]

Ah. Might the problem be that there's a 0 in the data? That 0 means you compute
–(a – 1) * ∞
which is big trouble for any value of a that's not 1. On the
unconstrained scale, there's actually no valid setting of a. But even
on the constrained scale, the prior assigns zero probability to a=1,
so you're guaranteed that the log-density of the data will be either
infinity or negative infinity. What happens if you replace the 0 with
0.001?

Matt

[Bob Carpenter]

After some discussion on the dev list, Matt correctly
guessed there was a zero in the data being modeled as
beta distributed, which I repeat here, with the zero in the
first line of data highlighted.

O <-
structure(c(0.7, 0.2, 0.55, 0.2, 0.5, 0, 0.42, 0.4, 0.3, 0.8, 0.74, 0.6, 0.72,
^
0.19, 0.58, 0.18, 0.48, 0.01, 0.45, 0.48, 0.4, 0.85, 0.75, 0.68, 0.69, 0.26,
0.54, 0.19, 0.55, 0.05, 0.44, 0.49, 0.35, 0.78, 0.71, 0.69, 0.68, 0.24, 0.52,
0.21, 0.51, 0.04, 0.46, 0.5, 0.38, 0.81, 0.7, 0.7, 0.71, 0.35, 0.51, 0.22, 0.52,
0.03, 0.43, 0.51, 0.39, 0.82, 0.72, 0.71),
.Dim = c(4, 3, 5))

I don't know there's much we can do to fix this issue. An actual 0 is
outside the proper support for the beta distribution.

- Bob

[Bob Carpenter]

We could pretty easily add a version of beta for logit-scaled parameters
directly,

logit_beta(y|a,b) = beta(inv_logit(y)|a,b) * |d/du inv_logit(u)|_y,

which is arithmetically stable for y after a bit of algebra,

log logit_beta(y|a,b) = alpha * log_sum_exp(0,y)
+ beta * log_sum_exp(0,-y)
+ lgamma(a + b) - lgamma(a) - lgamma(b).

But what should we do if a or b is close to 0?

We have a further issue when a and b in beta(a,b) are themselves
parameters, becuase what's really going on is beta(exp(a'),exp(b'))
where a' and b' the unconstrained parameters being sampled by HMC/NUTS.

If inv_logit(y) is used as a probability parameter in
bernoulli or binomial, we can deal with it directly using bernoulli_logit
and binomial_logit. So an explicit beta binomial model would be

logit_theta ~ logit_beta(a,b);
n ~ binomial_logit(N,logit_theta);

where inv_logit(logit_theta) has the same distribution as theta under
the standard parameterization,

theta ~ beta(a,b);
n ~ binomial(N,y);

The logit version should be more efficient and more arithmetically
stable because we don't have to do any transforms of the parameters.

- Bob

[Matt Hoffman]

What to do when a or b are very small is another question, one that's
relevant to either beta or beta_logit. Maybe there's something we
could do with Stirling's approximation to the (log of the) beta
function?
http://en.wikipedia.org/wiki/Beta_function

Matt