Beta function and integral value

Question

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

Answer 1

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):
    return integrate.quad(lambda t: t**(a-1)*(1-t)**(b-1),0,1)[0]

def beta2(a,b):
    return 2*integrate.quad(lambda theta: np.sin(theta)**(2*a-1)*np.cos(theta)**(2*b-1),0,np.pi/2)[0]

def beta3(a,b):
    return integrate.quad(lambda q: q**(a-1)/(1+q)**(a+b),0,np.infty)[0]

def func1(a,b):
    return 0.5*integrate.quad(lambda p: abs(p**(a-1)*(1-p)**(b-1)/(beta1(a,b))-1),0,1)[0]

def func2(a,b):
    return 0.5*integrate.quad(lambda p: abs(p**(a-1)*(1-p)**(b-1)/(beta2(a,b))-1),0,1)[0]

def func3(a,b):
    return 0.5*integrate.quad(lambda p: abs(p**(a-1)*(1-p)**(b-1)/(beta3(a,b))-1),0,1)[0]

def func4(a,b):
    return 0.5*integrate.quad(lambda p: abs(p**(a-1)*(1-p)**(b-1)/(special.beta(a,b))-1),0,1)[0]

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):
    return 0.5*integrate.quad(lambda p: abs(p**(a-1)*(1-p)**(b-1)/(betastirling(a,b))-1),0,1)[0]

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
The value using beta_2: 0.9204613256580434
The value using scipy: 0.9092842735473163