09-05-2021, 05:59 AM
Hello, below is a program to calculate the incomplete regularized beta function. The program is written with PYTHON (this program is available on the Web). When I launch the program, I only get a zero on the screen. Where am I doing wrong?
Thanks for your cooperation, Roberto.
Thanks for your cooperation, Roberto.
Code:
#PYTHON EXPORT nome
from math import *
def contfractbeta(a,b,x, ITMAX = 200):
""" contfractbeta() evaluates the continued fraction form of the incomplete Beta function; incompbeta().
(Code translated from: Numerical Recipes in C.)"""
EPS = 3.0e-7
bm = az = am = 1.0
qab = a+b
qap = a+1.0
qam = a-1.0
bz = 1.0-qab*x/qap
for i in range(ITMAX+1):
em = float(i+1)
tem = em + em
d = em*(b-em)*x/((qam+tem)*(a+tem))
ap = az + d*am
bp = bz+d*bm
d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem))
app = ap+d*az
bpp = bp+d*bz
aold = az
am = ap/bpp
bm = bp/bpp
az = app/bpp
bz = 1.0
if (abs(az-aold)<(EPS*abs(az))):
return az
def incompbeta(a, b, x):
''' incompbeta(a,b,x) evaluates incomplete beta function, here a, b > 0 and 0 <= x <= 1. This function requires contfractbeta(a,b,x, ITMAX = 200)
(Code translated from: Numerical Recipes in C.)'''
if (x == 0):
return 0;
elif (x == 1):
return 1;
else:
lbeta = lgamma(a+b) - lgamma(a) - lgamma(b) + a * log(x) + b * log(1-x)
if (x < (a+1) / (a+b+2)):
return exp(lbeta) * contfractbeta(a, b, x) / a;
else:
return 1 - exp(lbeta) * contfractbeta(b, a, 1-x) / b;
#end
EXPORT BetaRegInc(a,b,x)
BEGIN
PYTHON(nome,a,b,x)
END;