# Beta function and integral value

I have two values $$a$$ and $$b$$ where $$a \ge 0$$ and $$b \ge 0$$ and I have to calculate the formula below.

$$\frac{1}{2}\int_0^1\text{abs}\left[\left( \frac{p_i^{(a - 1)} \times (1 - p_i)^{(b - 1)}}{\beta(a,b)} \right)- 1\right] dp_i$$ I have to calculate above value. I am facing below error while doing that.

...

My equation is implemented as

import scipy.special as sc
import scipy.integrate as spi
integrand_np = lambda p_i:   abs(((p_i ** (a - 1)) * ((1 - p_i) ** (b - 1)) / (sc.beta(a,b)) - 1) / 2

result1, error = spi.quad(integrand_np, a, b)


:35: IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated. result1, error = spi.quad(integrand_np, a, b) ...

I am getting an error. I tried to use betaln by changing my equation. Still I am facing below issue

IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated.

Sample a values = [1400, 440, 799] Sample b values = [700,560, 1400]

To prevent error is there a way to use scipy.special.betaln. I changed above formula as below to use betaln

$$\frac{1}{2}\int_0^1 ( \frac{p_i^{(a - 1)} \times (1 - p_i)^{(b - 1)})}{\beta(a,b)} ) dp_i - \frac{1}{2}\int_0^1 1 dp_i$$

$$\frac{1}{2*beta(a,b)}\int_0^1 ( p_i^{(a - 1)} \times (1 - p_i)^{(b - 1)}) ) dp_i - \frac{1}{2}\int_0^1 1 dp_i$$

$$integralvalue = \int_0^1 ( p_i^{(a - 1)} \times (1 - p_i)^{(b - 1)}) ) dp_i$$

Then I am applying log to use betaln function. But still it is not working. ...

import scipy.special as sc
import scipy.integrate as spi
import numpy as np

integrand = lambda p_i: abs(((p_i ** (p_count - 1)) * ((1 - p_i) ** (n_count - 1))))

integral_second_term = lambda p_i: 1/2
result1, error = spi.quad(integrand, a, b)

result2, error1 = spi.quad(integral_second_term, a, b)
final_log = np.log(result1) - sc.betaln(p_count, n_count) - np.log(2)
final_output = abs((np.exp(final_log)) - result2)


...

It is not giving correct values

• Please use mark-up to show your code in a reasonably readable format. It would also help if you showed us what the error you get actually is! Commented Dec 15, 2022 at 15:58
• Yes $p_{i}$ in [0,1] range Commented Dec 20, 2022 at 20:00
• For $a$ and $b$ of the magnitude indicated in the question, $p_{i}^{(a-1)}(1-p_{i})^{(b-1)}$ underflows to zero when the computation is performed in IEEE-754 double precision. This is presumably why SciPy complains. It would probably be best to tackle this at the algorithm level: For what purpose is this integral being evaluated, and what alternative computations could be used instead? Commented Dec 20, 2022 at 22:39
• For what it's worth, $\frac{p_{i}^{(a-1)} (1 - p_{i})^{(b-1)}}{\mathrm{B}(a,b)} = \frac{\exp \left((a-1) \log (p) + (b-1) \log (1 - p) + \left(\frac{1}{2} - a\right) \log (a) + \left(\frac{1}{2} - b\right) \log (b) + \left(a+b-\frac{1}{2}\right) \log (a + b)\right)}{g \sqrt {2 \pi}}$, where $g=\frac{\Gamma^{\star}(a)\Gamma^{\star}(b)}{\Gamma^{\star}(a+b)}$, and $\Gamma^{\star}$ is the regulated gamma function (gammastar) introduced by Temme. For the magnitude of $a$, $b$ indicated in the question, this likewise underflows to zero. Commented Dec 20, 2022 at 22:55
• @BhavanaReddy Isn't what is needed here simply the CDF of the beta distribution, available in SciPy as scipy.special.btdtr(a, b, x) (here with $x=1$)? Commented Dec 21, 2022 at 3:11

I've tried to implement your equation \begin{align} f(a,b) = \frac{1}{2}\int_0^1 \left|\frac{p^{a-1}(1-p)^{b-1}}{\beta(a,b)}-1\right| \, \text{d}p, \end{align} using four different methods of computing the beta function (Euler's beta integral). These are given by \begin{align} \beta_1(a,b) &= \int_0^1 t^{a-1}(1-t)^{b-1} \, \text{d}t \\ \beta_2(a,b) &= 2\int_0^{\pi/2} \sin^{2a-1}\theta\cos^{2b-1}\theta \, \text{d}\theta \\ \beta_3(a,b) & = \int_0^\infty \frac{t^{a-1}}{(1+t)^{a+b}} \, \text{d}t, \end{align} and $$\beta_4$$ is scipy.special.beta(a,b) (Scipy link). These are implemented as

import numpy as np
from scipy import integrate
from scipy import special

def beta1(a,b):

def beta2(a,b):

def beta3(a,b):

def func1(a,b):

def func2(a,b):

def func3(a,b):

def func4(a,b):


This seems to work for reasonable values for $a$ and $b$. As an example for $a=15, b=15$ the results are

The value using beta_1: 0.605917291343624
The value using beta_2: 0.6059172913434613
The value using beta_3: 0.6059172913435065
The value using scipy: 0.605917291343624


For your sets of parameters I only managed to get results for $a=440, b=560$. For large $a$ and $b$ one can apply Stirling's approximation to the beta function (Wiki link) \begin{align} \beta(a,b) \simeq \sqrt{2\pi} \frac{a^{a-1/2}b^{b-1/2}}{(a+b)^{a+b-1/2}}, \end{align} implemented as

def betastirling(a,b):
return np.sqrt(2*np.pi)*a**(a-0.5)*b**(b-0.5)/(a+b)**(a+b-0.5)

def funcstirling(a,b):

Which again works for reasonable values of $a$ and $b$. Here are the results for $a=440, b=560$
The value using beta_1: 0.9103593870524365