1
$\begingroup$

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

$\endgroup$
7
  • $\begingroup$ 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! $\endgroup$ Dec 15, 2022 at 15:58
  • $\begingroup$ Yes $p_{i} $ in [0,1] range $\endgroup$ Dec 20, 2022 at 20:00
  • $\begingroup$ 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? $\endgroup$
    – njuffa
    Dec 20, 2022 at 22:39
  • $\begingroup$ 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. $\endgroup$
    – njuffa
    Dec 20, 2022 at 22:55
  • $\begingroup$ @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$)? $\endgroup$
    – njuffa
    Dec 21, 2022 at 3:11

1 Answer 1

3
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ Thanks for equations. 440 and 560 special.beta is only working. $\endgroup$ Dec 21, 2022 at 18:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.