Solved – Beta distribution fitting in Scipy

beta distributionpythonscipy

According to Wikipedia the beta probability distribution has two shape parameters: $\alpha$ and $\beta$.

When I call scipy.stats.beta.fit(x) in Python, where x is a bunch of numbers in the range $[0,1]$, 4 values are returned. This strikes me as odd.

After googling I found one of the return values must be 'location', since the third variable is 0 if I call scipy.stats.beta.fit(x, floc=0).

Does anyone know what the fourth variable is, and if the first two are $\alpha$ and $\beta$?

Best Answer

Despite an apparent lack of documentation on the output of beta.fit, it does output in the following order:

$\alpha$, $\beta$, loc (lower limit), scale (upper limit - lower limit)

Related Solutions

Solved – How to sample from truncated distributions using scipy

Perhaps the simplest and most generic method is inverse transform sampling. All of the distributions mentioned in OP are absolutely continuous, so it's very straightforward. Morevoer, for these distributions, the key components are all implemented in scipy already, so we only need to do algebra.

Suppose $X$ is distributed according to some CDF $G$. The truncated distribution $F$ is how $x$ is distributed given that it's restricted to the interval $[a,b]$. This is just rescaling and shifting the CDF $G$, so we have $$ F(y) = \frac{G(y) - G(a)}{G(b)-G(a)}. $$

Inverse transform sampling observes that for some continuous random variable, we can sample from a CDF $F$ using a uniform distribution. This is because $F(F^{-1}(u))=u$.

Putting it all together,

$$ F^{-1}(u)=G^{-1}\left[u(G(b) - G(a)) +G(a) \right] $$

scipy implements CDFs $G$ and quantile functions $G^{-1}$ for all of the distributions you list, so the rest is just writing some (simple!) code:

draw a value of $u$ from a Uniform(0,1) distribution
compute $F^{-1}(u)$

There are lots of other algorithms to solve this problem (rejection sampling, ziggurat algorithm and many more), at various trade-offs of speed, precision and complexity of implementation. I think this one sits at a nice compromise.

If this seems too complicated, then you could do something even more naive and just reject any values that are outside of your desired bounds. This is inefficient, especially if $F(b) - F(a)$ the probability in the interval $[a,b]$ is small.

Bayesian Multilevel Models – Inconsistent Posterior Estimates in Beta-Binomial Likelihood vs. Binomial

I posted a similar question at the PyMC3 discourse. The answer was to do proper multilevel modeling/pooling. Essentially, I should have one $\theta$ per observation and each one should be "communicate" with one another through the common $\alpha$ and $\beta$ parameters. The results were identical to having used the Beta-Binomial likelihood.

It seems that trying to infer a single $\theta$, it wasn't able to vary as knowledge of $\alpha, \beta$ determined it, exactly (and conversely.) But by turning $\theta$ into a vector, rather than a scalar, I got good inference on $\alpha, \beta$.

with pm.Model() as m1:
  alpha = pm.HalfNormal('alpha',sigma=1)
  beta = pm.HalfNormal('beta',sigma=1)
  theta = pm.Beta('theta',alpha, beta, shape=len(W))
  lik = pm.Binomial('likelihood',n=N,p=theta,observed=W)
  trace1 = pm.sample(chains=4)

Best Answer

Related Solutions

Solved – How to sample from truncated distributions using scipy

Bayesian Multilevel Models – Inconsistent Posterior Estimates in Beta-Binomial Likelihood vs. Binomial

Related Question