where B(α, β) is the beta
function, implemented in this library as beta.
Division by the beta function ensures that the pdf is normalized to the
range zero to unity.

The following graph illustrates examples of the pdf for various values
of the shape parameters. Note the α = β = 2 (blue line) is dome-shaped, and
might be approximated by a symmetrical triangular distribution.

If α = β = 1, then it is a ​
uniform
distribution, equal to unity in the entire interval x = 0 to
1. If α ​ and β ​ are < 1, then the pdf is U-shaped. If α != β, then the shape
is asymmetric and could be approximated by a triangle whose apex is away
from the centre (where x = half).

Requires alpha,beta > 0,otherwise domain_error
is called. Note that technically the beta distribution is defined for
alpha,beta >= 0, but it's not clear whether any program can actually
make use of that latitude or how many of the non-member functions can
be usefully defined in that case. Therefore for now, we regard it as
an error if alpha or beta is zero.

For example:

beta_distribution<>mybeta(2,5);

Constructs a the beta distribution with alpha=2 and beta=5 (shown in
yellow in the graph above).

One estimates either α ​ or β ​
from presumed-known mean and variance.

The other pair estimates either α ​ or β ​ from the cdf and x.

It is also possible to estimate α ​ and β ​ from 'known' mode & quantile.
For example, calculators are provided by the Pooled
Prevalence Calculator and Beta
Buster but this is not yet implemented here.

With integer values of α ​ and β ​ the distribution B(i, j) is that of the j-th
highest of a sample of i + j + 1 independent random variables uniformly
distributed between 0 and 1. The cumulative probability from 0 to x is
thus the probability that the j-th highest value is less than x. Or it
is the probability that that at least i of the random variables are less
than x, a probability given by summing over the Binomial
Distribution with its p parameter set to x.