In R there is the beta.select() function, which finds the shape parameters of a beta density that matches the knowledge of two quantiles of the distribution.

Here is an example from the documentation:

quantile1=list(p=.5,x=0.25)
quantile2=list(p=.9,x=0.45)
beta.select(quantile1,quantile2)

However, to do the same in Python, I haven't found any similar method in the beta documentation of scipy.

• Is there any alternative to the scipy library to perform this task?
• What is the equivalent of the R's quantile generation with list() in Python? I know how to compute the quantile of a given array with np.quantile(), but it expect an array as input, not the probability (e.g p=.5).

CodePudding user response：

There is no equivalent function in SciPy. Here's an implementation based on scipy.special.betainc and scipy.optimize.fsolve:

from scipy.special import betainc
from scipy.optimize import fsolve


def _beta_select_equation(params, p1, x1, p2, x2):
    return betainc(*params, [x1, x2]) - [p1, p2]


def beta_select(p1, x1, p2, x2):
    params, info, status, mesg = fsolve(_beta_select_equation, [1, 1],
                                        args=(p1, x1, p2, x2), xtol=1e-12,
                                        full_output=True)
    if status != 1:
        raise RuntimeError(f'fsolve failed: {mesg}')
    return params

For example,

In [21]: beta_select(p1=0.5, x1=0.25, p2=0.9, x2=0.45)
Out[21]: array([2.66897386, 7.36479059])