Regularisierte Betafunktion / Rekursion

[Bob Billsworth]

Hallo Zusammen

Ich möchte diese Repräsentation der regularisierten Betafunktion implementieren. Jede der von mir durchgegangenen Varianten -selbst im Iterationsverfahren- liefert für z = 0.5, a = 2.5 und b = 5 den Wert 0.2772238763345327. Das richtige Ergebnis lautet jedoch 0.83580505... .

Dies ist mein aktueller Stand:

Code: Alles auswählen

import math as ma

def betafunktion(x, y):
    return (ma.gamma(x)*ma.gamma(y))/ma.gamma(x+y)

def regbeta(z,a,b,c1=0,k=0):
    if c1 == 0:
        const = ((z**a)*((1-z)**b))/(a*betafunktion(a,b))
        var = 1+regbeta(z,a,b,c1+1,k+1)
        #print(const, var)
        return const/var
    elif c1 == 1:
        c1 = 2
        zaehler = -1*(((a+k)*(a+b+k)*z)/((a+2*k)*(a+2*k+1)))
        nenner = 1+regbeta(z,a,b,c1,k+1)
        #print("U: ", zaehler, nenner, zaehler/nenner, k)
        return zaehler/nenner
    else:
        if k == 200:
            return 10**-10
        else:
            c1 = 1
            zaehler = (k*(b-k)*z)/((a+2*k-1)*(a+2*k))
            nenner = 1+regbeta(z,a,b,c1,k+1)
            #print("G: ", zaehler, nenner, zaehler/nenner, k)
            return zaehler/nenner

Für Eure Hilfe danke ich Euch im Voraus.

[Bob Billsworth]

Ich habe mehrere Lösungen gefunden.

Die oben genannte Repräsentation der regularisierten Betafunktion lässt sich als -schnellere- iterative Variante implementieren. Der C-Code wird in "Numerical Recipes in C" beschrieben (siehe hier).

Präziser und schneller ist jedoch die Formel von Didonato und Morris. Ich habe diese iterativ gemäß dem von Thompson und Barnett ("Coulomb and Bessel Functions of Complex Arguments and Order", Journal of Computational Physics: Vol. 64, Nr.2, Juni 1986, Seite 508) modifizierten Lentz-Algorithmus sowie rekursiv implementiert.

Die etwas ungelenke rekursive Variante:

Code: Alles auswählen

import math as ma

def betafunktion(x, y):
    return (ma.gamma(x)*ma.gamma(y))/ma.gamma(x+y)

def regbetaDMREK(x,a,b,n=0):
    if n == 0:
        const = ((x**a)*((1-x)**b))/betafunktion(a,b)
        alpha0 = 1
        beta0 = (a/(a+1))*((a-(a+b)*x)+1)
        return const*(alpha0/(beta0+regbetaDMREK(x,a,b,n+1)))
    elif n > 0 and n < 600:
        betan = n+((n*(b-n)*x)/(a+2*n-1))+((a+n)/(a+2*n+1))*((a-(a+b)*x)+1+n*(1+(1-x)))
        alphan = (((a+n-1)*(a+b+n-1))/((a+2*n-1)**2))*(n*(b-n)*x**2)
        return alphan/(betan+regbetaDMREK(x,a,b,n+1))
    else:
        return 10**-15

Die sehr präzise iterative Variante:

Code: Alles auswählen

import math as ma

def betafunktion(x, y):
    return (ma.gamma(x)*ma.gamma(y))/ma.gamma(x+y)

def regbetaDMMLA(x,a,b):
    const = ((x**a)*((1-x)**b))/betafunktion(a,b)
    ergebnis = (a/(a+1))*((a-(a+b)*x)+1)
    D = 0
    C = ergebnis
    for n in range(1,2001):
        betan = n+((n*(b-n)*x)/(a+2*n-1))+((a+n)/(a+2*n+1))*((a-(a+b)*x)+1+n*(1+(1-x)))
        alphan = (((a+n-1)*(a+b+n-1))/((a+2*n-1)**2))*(n*(b-n)*x**2)
        D = betan + alphan * D
        C = betan + alphan / C
        D = 1/D
        delta = D * C
        ergebnis *= delta
        if ma.fabs(delta-1) < 10**-16:
            break
    return const/ergebnis

Wem an einem ganz stabilen Code gelegen ist, der sollte noch für "ergebnis", "D" und "C" jeweils die Überprüfung auf Null mit entsprechender Zuweisung eines sehr kleinen Werts einfügen.