# Beta distribution and floating point precision

**URL:** <https://discourse.mc-stan.org/t/beta-distribution-and-floating-point-precision/26485>\
**Category:** Developers\
**Created:** [February 21, 2022, 7:05am UTC](https://discourse.mc-stan.org/t/beta-distribution-and-floating-point-precision/26485 "2022-02-21T07:05:36Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![saudiwin](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/saudiwin/32/1137_2.png) [@saudiwin](https://discourse.mc-stan.org/u/saudiwin)\
**Post date:** [February 21, 2022, 7:05am UTC](https://discourse.mc-stan.org/t/beta-distribution-and-floating-point-precision/26485/1 "2022-02-21T07:05:36Z")

</div>

Hi all -

I’ve been working a lot with beta regression lately, and I am having some issues using probabilities that are censored by floating point precision, i.e. very close to 0/1. It strikes me that the best approach to dealing with this problem is to have a beta distribution on log probabilities, rather than the (0,1) interval. Any kind of rounding can introduce some pretty big biases as the extreme values will have an outsize effect on estimates (see my paper [here](https://osf.io/preprints/socarxiv/2sx6y/)). In principle it doesn’t seem to hard to have a Beta distribution on \text{log } x rather than x \in (0,1), but I can’t find any examples of this in the literature. If anyone knows of one, or can see a way of straightforwardly calculating the beta distribution on \text{log } x, I’d appreciate any tips/pointers.

---

<div class="post-metadata">

**Author:** ![saudiwin](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/saudiwin/32/1137_2.png) [@saudiwin](https://discourse.mc-stan.org/u/saudiwin)\
**Post date:** [February 21, 2022, 7:28am UTC](https://discourse.mc-stan.org/t/beta-distribution-and-floating-point-precision/26485/2 "2022-02-21T07:28:44Z")

</div>

This new paper seems like one plausible option… a “beta prime” distribution on \frac{x}{1-x} that includes a regression version to predict the mean of \frac{x}{1-x}.

> **[A new regression model for positive random variables with skewed and long tail...](https://link.springer.com/article/10.1007/s40300-021-00203-y)**
>
> In this paper, we propose a regression model where the response variable is beta prime distributed using a new parameterization of this distribution that is indexed by mean and precision parameters. The proposed regression model is useful for...

---

<div class="post-metadata">

**Author:** ![saudiwin](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/saudiwin/32/1137_2.png) [@saudiwin](https://discourse.mc-stan.org/u/saudiwin)\
**Post date:** [February 21, 2022, 10:31am UTC](https://discourse.mc-stan.org/t/beta-distribution-and-floating-point-precision/26485/3 "2022-02-21T10:31:23Z")

</div>

OK third follow-up answering my own question :O.

Based on the paper above, we can also model a beta regression using the following log-likelihood function for a variable x \in (0,1), where y = \frac{x}{1-x}, \mu \> 0, \phi \>0.

```nohighlight
    real beta_prime_lpdf(real y, real mu, real phi) {

    // adapted from Bourguignon et al. (2021)
    // doi: https://doi.org/10.1007/s40300-021-00203-y

        return((mu*(1 + phi) - 1)*log(y) - (mu*(1 + phi) + phi + 2)*log1p(y) -
          lgamma(mu*(1 + phi)) - lgamma(phi + 2) +
          lgamma(mu*(1+ phi) + phi + 2));

    }

```

It is then straightforward to translate model predictions back to probabilities by using the inverse logit function as it’s essentially modeling the log-odds. The author of the paper has a helpful R package [BPmodel](https://cran.r-project.org/web/packages/BPmodel/index.html) that implements distribution functions.

I’m putting this into `ordbetareg`, but I think it could be of wider use in the Stan community given the issues with overflow/underflow in the `beta_proportion_lpdf` function.

---

<div class="post-metadata">

**Author:** ![rok\_cesnovar](https://avatars.discourse-cdn.com/v4/letter/r/7bcc69/32.png) [@rok\_cesnovar](https://discourse.mc-stan.org/u/rok_cesnovar)\
**Post date:** [February 21, 2022, 11:13am UTC](https://discourse.mc-stan.org/t/beta-distribution-and-floating-point-precision/26485/4 "2022-02-21T11:13:19Z")

</div>

Thanks @saudiwin!

I am guessing putting it into @spinkney’s [GitHub - spinkney/helpful\_stan\_functions](https://github.com/spinkney/helpful_stan_functions) repository would be great for starters as well (before we embark on adding it in Stan Math). What do you think Sean?

---

<div class="post-metadata">

**Author:** ![saudiwin](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/saudiwin/32/1137_2.png) [@saudiwin](https://discourse.mc-stan.org/u/saudiwin)\
**Post date:** [February 21, 2022, 1:41pm UTC](https://discourse.mc-stan.org/t/beta-distribution-and-floating-point-precision/26485/5 "2022-02-21T13:41:19Z")

</div>

Sounds good.

The one caveat with the above approach is that it can’t estimate all the possible shapes of the Beta distribution. The \alpha and \beta parameters are restricted to be at least +1 and +2 respectively. Apparently that’s because the mean and variance of the beta prime distribution aren’t defined otherwise.

It does work pretty well, though, given that caveat.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [March 9, 2022, 9:55pm UTC](https://discourse.mc-stan.org/t/beta-distribution-and-floating-point-precision/26485/6 "2022-03-09T21:55:19Z")

</div>

> [@saudiwin](#):
>
> In principle it doesn’t seem to hard to have a Beta distribution on \text{log } x log x\text{log } x rather than x \in (0,1) x∈(0,1)x \in (0,1)

Given \theta \sim \textrm{beta}(\alpha, \beta), then we can let \phi = \log \theta or \psi = \log(\theta / (1 - \theta)) and work out p(\phi \mid \alpha, \beta) and p(\psi \mid \alpha, \beta) in the usual way. So let’s consider \phi, for which the change-of-variables formula yields

\begin{eqnarray\*} \log p\_{\Phi}(\phi) & = & \log \left( \textrm{beta}(\exp(\phi) \mid \alpha, \beta) \cdot \left| \frac{d}{d\phi} \exp(\phi) \right| \right) \\[4pt] & = & \log \left( \exp(\phi)^{\alpha - 1} (1 - \exp(\phi))^{\beta - 1} \cdot \exp(\phi) \right) + \textrm{const}. \\[4pt] & = & (\alpha - 1) \cdot \phi + (\beta - 1) \cdot \textrm{log1m\_exp}(\phi) + \phi + \textrm{const}. \end{eqnarray\*}

If \alpha and \beta are unknowns, then you need to evaluate the \Gamma(\alpha + 1), \Gamma(\beta + 1) and \Gamma(\alpha + \beta + 1), which can also be problematic.

The original has

\log \textrm{beta}(\theta \mid \alpha, \beta) = (\alpha - 1) \cdot \theta + (\beta - 1) \cdot (1 - \theta) + \textrm{const}.

and the main risk if you don’t have the normalizer is that 1 - \theta rounds to 1.

Now if you’re using `beta_lpdf` rather than `beta_lupdf` or the sampling statement, you can also run into problems evaluating the beta function normalizer, which I’ve just rendered as const here.