Exakte Berechnung der unvollständigen Betafunktion Βₓ(a,b) und Iₓ(a,b)
mit ganzzahligen Koeffizienten und rationalem x

Die englische Wikipedia schreibt (ohne weitere Begründung):

For positive integer a and b, the incomplete beta function will be a polynomial of degree a + b - 1 with rational coefficients.

Ins Deutsche übersetzt:

⟪ Β_x(a,b) = sum_(i=0)^(a+b-1) k_i x^i " mit " k_i in QQ ⟫

Das prüfen wir doch einmal nach.

1. Umformung der Integral-Darstellung von Βₓ(a,b) mit ganzzahligen a und b zu einem Polynom

Schauen wir uns die Definition von Βₓ(a,b) an:

(1) ⟪ Β_x(a,b) = int_0^x t^(a-1) (1-t)^(b-1) dt ⟫

Um die Potenz von 1 - t zu beseitigen, nutzen wir die erweiterte binomische Formel:

(2) ⟪ (x + y)^m = sum_(i=0)^m ((m),(i)) \ x^(m-i) \ y^i ⟫

Und wenden sie mit x := 1, y := -t und m := b-1 an:

(3) ⟪ Β_x(a,b) = int_0^x t^(a-1) sum_(i=0)^(b-1) ((b-1),(i)) \ 1^(b-1-i) \ (-t)^i \ dt ⟫

⟪ Β_x(a,b) = int_0^x t^(a-1) sum_(i=0)^(b-1) ((b-1),(i)) \ (-t)^i \ dt ⟫

⟪ Β_x(a,b) = int_0^x t^(a-1) sum_(i=0)^(b-1) ((b-1),(i)) \ (-1)^i \ t^i \ dt ⟫

Wir ziehen das t ^a - 1 in die Summe und erhalten so eine Gleichung, in der nur noch ein Faktor von t abhängt:

(4) ⟪ Β_x(a,b) = int_0^x sum_(i=0)^(b-1) (-1)^i \ ((b-1),(i)) \ t^(a-1+i) \ dt ⟫

Nach Vertauschen von Integral und Summation können wir integrieren.

(5) ⟪ Β_x(a,b) = sum_(i=0)^(b-1) (-1)^i \ ((b-1),(i)) \ int_0^x t^(a-1+i) \ dt ⟫

Mit ∫xⁿ = xⁿ⁺¹ / (n+1) erhalten wir:

(6) ⟪ Β_x(a,b) = sum_(i=0)^(b-1) (-1)^i \ ((b-1),(i)) \ x^(a+i) / (a+i) ⟫

Substitution von i durch i-a ergibt:

⟪ Β_x(a,b) = sum_(i=a)^(a+b-1) (-1)^(i-a) \ ((b-1),(i-a)) \ 1/i \ x^i ⟫

(7) ⟪ Β_x(a,b) = sum_(i=0)^(a+b-1) k_i \ x^i " mit " k_i = { ( (-1)^(i-a) \ ((b-1),(i-a)) \ 1/i, if i ge a ), ( 0 , if i lt a ) :} ⟫

Et voilà, da ist die Darstellung als Polynom mit rationalen Koeffizienten; und die höchste Potenz von x ist offensichtlich a+b-1. ☺

2. Berechnung für x = 1

Für x = 1 vereinfacht sich die Formel (6) ein wenig:

(8) ⟪ Β_1(a,b) = sum_(i=0)^(b-1) (-1)^i ((b-1),(i)) 1/(a+i) ⟫

Schauen wir uns die Funktionswerte für kleine Werte von a und b an:

Die Funktion ist wie zu erwarten symmetrisch. Wir können so erweitern, dass der Nenner zu (a-b-1)! wird:

Die Zähler jeder Zeile betragen jeweils ein Vielfaches der Zähler der vorausgehenden Zeile, und die Zähler der Zeile a=1 haben jeweils den Wert (b-1)!. Versuchen wir also, im Zähler (b-1)! auszuklammern:

Offensichtlich ist der rot markierte erste Term jeweils gleich (a-1)!:

Die Funktionswerte für kleine a und b legen nahe, dass gilt:

(9) ⟪ Β_1(a,b) ≟ ((a-1)! (b-1)!) / ((a+b-1)!) ⟫

3. Beweis von Β₁(a,b) = (a-1)! * (b-1)! / (a+b-1)!

Wegen der Summe in der Darstellung von Β₁(a,b) bietet sich ein Induktionsbeweis über b an.

Wir verankern bei b = 1:

(10) ⟪ Β_1(a,[1]) = sum_(i=0)^([1]-1) (-1)^i (([1]-1),(i)) 1/(a+i) = (-1)^0 ((0),(0)) 1/a = 1/a = ((a-1)! ([1]-1)!) / ((a+[1]-1)!) ✓ ⟫

Für den Induktionsschritt nutzen wir eine Eigenschaft des Binomialkoeffizienten:

⟪ ((n+1),(k+1)) = ((n),(k+1)) + ((n),(k)) ⟫

Dabei erlauben wir k = -1, um ein unübersichtliches Herumgewürge mit Indizes zu vermeiden:

⟪ ((n),(-1)) := 0 ⟫

⟪ 1 = ((n+1),(0)) = ((n),(0)) + ((n),(-1)) = 1 + 0 = 1 ⟫

Für b + 1 erhalten wir:

(11) ⟪ Β_1(a,b+1) = sum_(i=0)^(b) (-1)^i ((b),(i)) 1/(a+i) ⟫

(12) ⟪ Β_1(a,b+1) = sum_(i=0)^(b) (-1)^i ((b-1),(i)) 1/(a+i) + sum_(i=0)^(b) (-1)^i ((b-1),(i-1)) 1/(a+i) ⟫

In der linken Summe können wir den letzten Summanden und in der rechten Summe den ersten Summanden streichen, weil der jeweilige Binomialkoeffizient 0 ist:

⟪ Β_1(a,b+1) = sum_(i=0)^(b-1) (-1)^i ((b-1),(i)) 1/(a+i) + sum_(i=1)^(b) (-1)^i ((b-1),(i-1)) 1/(a+i) ⟫

In der rechten Summe substituieren wir i durch i + 1:

⟪ Β_1(a,b+1) = sum_(i=0)^(b-1) (-1)^i ((b-1),(i)) 1/(a+i) + sum_(i=0)^(b-1) (-1)^(i+1) ((b-1),(i)) 1/(a+(i+1)) ⟫

Wir ziehen ein -1 aus der rechten Summe:

⟪ Β_1(a,b+1) = sum_(i=0)^(b-1) (-1)^i ((b-1),(i)) 1/(a+i) - sum_(i=0)^(b-1) (-1)^i ((b-1),(i)) 1/((a+1)+i) ⟫

Und haben damit eine Rekursionsformel:

(13) ⟪ Β_1(a,b+1) = Β_1(a,b) - Β_1(a+1,b) ⟫

Jetzt müssen wir nur noch einsetzen:

⟪ Β_1(a,[b+1]) = ((a-1)! \ (b-1)!) / ((a+b-1)!) - (([a+1]-1)! \ (b-1)!) / (([a+1]+b-1)!) ⟫

⟪ Β_1(a,[b+1]) = ((a-1)! \ (b-1)! \ (a+b)) / ((a+b-1)! \ (a+b)) - (a! \ (b-1)!) / ((a+b)!) ⟫

⟪ Β_1(a,[b+1]) = ((a-1)! \ (b-1)! \ (a+b) - a \ (a-1)! \ (b-1)!) / ((a+b)!) ⟫

⟪ Β_1(a,[b+1]) = ((a-1)! \ (b-1)! \ (a+b - a )) / ((a+b)!) = ((a-1)! \ b! ) / ((a+b)!) ⟫

⟪ Β_1(a,[b+1]) = ((a-1)! \ ([b+1]-1)! ) / ((a+[b+1]-1)!) ✓ ⟫

Zusammen mit (10) haben wir (9) bewiesen. ☺

4. Effiziente Berechnung

Wenn wir für mehrere x_j gleichzeitig rechnen, brauchen wir C(i) und p(i) nur einmal zu bestimmen:

(14) ⟪ Β_(x_j)(a,b) = sum_(i=0)^(b-1) (-1)^i \ ubrace( ubrace( ((b-1),(i)))_(C(i)) \ 1/(a+i) )_(p(i)) \ ubrace( x_j^(a+i))_(X_j(i)) ⟫

Die Faktoren C(i) und X_j(i) können inkrementell berechnet werden:

Um gemeinsame Ausdrücke nur einmal zu berechnen, führen wir eine Hilfsvariable ein:

(16) ⟪ p(i) = C(i) // (a+i) ⟫

Und summieren folgende Terme auf:

(17) ⟪ s_j(i) = p(i) \ X_j(i) ⟫

Das ergibt folgenden Algorithmus:

function beta( x[], a, b ){

	// Berechne C(0) und Xj(0)
	C  := 1
	Xj := pow(xj, a)
	i  := 0

	// Berechne p(0) und initialisiere Βj(0)
	p  := C / (a+i)
	Βj := p * Xj

	while( i < b-1 ){

		// Berechne C(i+1) und Xj(i+1)
		C  := C * (b-1-i) / (i+1)
		Xj := Xj * xj
		i  := i + 1

		// Berechne p(i) und aktualisiere Βj(i)
		p  := C / (a+i)
		if odd(i)
			Βj := Βj - p * Xj
		else
			Βj := Βj + p * Xj
		endif
	}

	return Β[]
}

5. Online-Ausführung von Algorithmus 1

Bei Werten ab etwa 10 für a und b führt diese naive Implementierung mit double-Genauigkeit zu unakzeptablen bis absurden Rechenfehlern, verursacht durch extreme Unterschiede der Größenordnungen der Summanden.
Nicht nachmachen!

Weil die Werte von Β_x(a,b) rasch sehr klein werden, wird auch der normalisierte Wert I_x(a,b)=Β_x(a,b)/Β₁(a,b) angezeigt.

6. Exakte Berechnung für große Werte von a und b und rationale Werte von x

Für rationale Werte von x ist der Wert des Polynomes (6) als Summe von Produkten rationaler Zahlen wieder rational und damit exakt berechenbar.

A straightforward binomial expansion of the integrand and a subsequent term-by-term integration results in an alternating series in powers of x which cannot be used for large values of a. The eventual subtraction of consecutive terms of nearly equal absolute values causes a prohibitive loss in significant digits.

A. R. DiDonato and M. P. Jarnagin (1967)

Hehe, doch genau das machen wir jetzt.

Dabei vermeiden wir einen Signifikanz-Verlust, indem wir mit exakten rationalen und natürlichen Zahlen rechnen, und bekommen so ein exaktes Ergebnis. Wobei dieses im allgemeinen ein Bruch mit sehr großem Zähler und Nenner sein wird. 🤪

Die rationale Zahlen werden als Produkt von Primzahlpotenzen mit ganzzahligem, auch negativem Exponent dargestellt:

⟪ q = prod_i p_i^(e_i) " mit " p_i in bbb "P", e_i in ZZ⟫

Beispiele: 20 = 2² · 5¹, 24 = 2³ · 3¹, 100 = 2² · 5², 1/2 = 2^-1, 2/3 = 2¹ · 3^-1, 4/9 = 2² · 3^-2.

Die ganzzahlige Werte i, i+1, a+i und b-1-i des Algorithmus werden vor Verwendung durch Primfaktorzerlegung in die Darstellung als Produkt von Primzahlpotenzen konvertiert.

Multiplikationen und Divisionen von Zahlen in dieser Darstellung sind sehr schnell, wobei das Kürzen automatisch erfolgt, anders als bei der Darstellung als Paar von zwei beliebig genauen ganzen Zahlen. Das gilt auch für die Berechnung des kleinsten gemeinsamen Vielfaches oder des größten gemeinsamen Teiler.

Leider lassen sich Zahlen in dieser Darstellung nicht addieren oder subtrahieren, diese müssen zuerst in ein Darstellung als Polynom überführt werden, in der Addition und Subtraktionen trivial sind:

⟪ n = sum_j d_j K^j " mit " d_j in NN_0, 0 le d_j lt K⟫

Und in der Tat braucht diese Konversion, die implizit bei der Summation über S[i] * N im zweiten Teil des Algorithmus durchgeführt wird, den größten Teil der Laufzeit.

function beta( x, a, b ){

	// (1) Berechne Summanden als exakte rationale Zahlen S = Sz/Sn

	// Berechne C(0) und X(0)
	C := 1        // oder (a+b-1)!/(a-1)!/(b-1!) für Berechnung von Iₓ(a,b)
	X := pow(x, a)
	i := 0

	// Berechne und speichere Summanden
	S[0] := C / (a+i) * X

	while( i < b-1 ){

		// Berechne C(i+1) und Xj(i+1)
		C := C * (b-1-i) / (i+1)
		X := X * x
		i := i + 1

		// Berechne und speichere Summanden
		S[i] := C / (a+i) * X
	}

	// (2) Mache die Summanden ganzzahlig und addiere sie

	// bestimme kleinstes gemeinsames Vielfaches N der Nenner der Summanden
	N := kgvi( numerator(Si) )

	// multipliziere Summanden mit N, addiere die (ganzzahligen) Produkte
	Z := ∑i( S[i] * N )

	// betaₓ(a,b) ist exakt Z/N
	return ( Z, N )
}

Die Dezimaldarstellung des Bruch bekommt man durch Multiplikation des Zählers Z mit der passenden Zehnerpotenz, ganzzahliger Division des Produktes durch den Nenner N und Einfügen des Dezimalpunktes. Um einen gerundeten Wert zu bekommen, addiert man vor der Division die Hälfte des Nenners zum Zähler.

Den Wert der regularisierten Beta-Funktion I_x(a,b) bekommt man durch Division von Β_x(a,b) durch Β₁(a,b).

7. Online-Ausführung der exakten Berechnung

Benutzt ganze Zahlen und rationale Zahlen mit unbegrenzter Genauigkeit.

Eine Berechnung kann sehr lange dauern, deshalb wird wird sie von einem Webworker durchgeführt und es wird ein Fortschrittsbalken angezeigt.

Da immer zuerst das exakte Ergebnis bestimmt wird, hat die gewählte Zahl signifikanter Ziffern der Dezimaldarstellung praktisch keinen Einfluss auf die Rechenzeit.