Die Herleitung der Stirling-Approximation mit der Laplace-Methode

Einordnung des Artikels

• Ausgewählte Kapitel der Mathematik (für Programmierer, Informatiker, Ingenieure und Naturwissenschaftler)
  • Spezielle Kapitel der Analysis
    • Die Herleitung der Stirling-Approximation mit der Laplace-Methode

Als mathematische Voraussetzungen werden hier verwendet, aber nicht erläutert: Integralrechnung, insbesondere uneigentliche Integrale, die Substitutionsregel und Gauß-Integrale.

Letztere werden in Eigenschaften von Zufallsvariablen: Die Varianz und die Standardabweichung ausführlich behandelt.

Einführung

Möchte man die Fakultäten n! von großen natürlichen Zahl n berechnen, so ist man häufig nicht an dem konkreten Zahlenwert interessiert, sondern an dem asymptotischen Verhalten von n!. Damit ist gemeint, dass man eine kontinuierliche Funktion F(n) sucht (F soll hier an Fakultät erinnern), die für große n ein ähnliches Verhalten zeigt wie n!. Berechnet man n! für kleine Werte von n (etwa bis zur Größenordnung von n = 100), so ist nicht klar, welches asymptotische Verhalten die Funktion F(n) haben wird:

• Wächst sie mit einer Potenz von n? Wenn ja mit welcher?
• Wächst sie schneller als mit einer Potenz, zum Beispiel exponentiell?
• Oder wächst sie sogar schneller als eine Exponentialfunktion?

Die bekannteste Näherung für die Fakultät n! ist die Stirling-Approximation, die eine Funktion F(n) angibt, die man auch für n weit jenseits der Größenordnung n = 100 leicht berechnen kann. Es gibt zahlreiche Möglichkeiten die Stirling-Approximation F(n) der Fakultät zu berechnen. Hier wird die sogenannte Laplace-Methode gezeigt, die weitreichende Anwendungen besitzt. Sie wurde von Laplace entwickelt, um Näherungen für eine Laplace-Transformation zu berechnen. Die hier gezeigte Herleitung der Stirling-Approximation ist somit ein Spezialfall dieser allgemeineren Methode.

Die Laplace-Methode (im allgemeinenen Sinn) wird hier nicht ausführlich vorgestellt. Sie wird nur verwendet, um die Stirling-Approximation herzuleiten. Ebenso wird hier nicht diskutiert, in welchem Sinn die Näherungsfunktion F(n) die Fakultät n! approximiert. Dazu müssten die entsprechenden Konvergenzbegriffe eingeführt und erläutert werden, was hier nicht geschehen soll.

Übersicht über die Vorgehensweise

Die Herleitung der Funktion F(n), die asymptotisch mit der Fakultät n! übereinstimmt, geschieht in mehreren Schritten:

1. Die Fakultät wird mit Hilfe der Gamma-Funktion Γ(x) ausgedrückt. Die Gamma-Funktion ist nicht nur für natürliche n, sondern für reelle x definiert; sie verallgemeinert somit die Fakultät. Für natürliche Zahlen n gilt: Γ(n+1) = n!.
2. Die Gamma-Funktion Γ(x) ist definiert durch ein Integral, das vom Parameter x abhängt; der Integrand ist das Produkt einer Potenzfunktion und einer Exponentialfunktion.
3. Durch geeignete Substitutionen wird das Integral zur Berechnung von n! so umgeformt, dass nur noch eine Exponentialfunktion über y integriert wird, nämlich exp(-n·f(y)). Dabei ergibt sich die Funktion f(y) aus den genannten Substitutionen.
4. Auf dieses Integral wird die sogenannte Laplace-Methode angewendet; diese besteht aus den folgenden Schritten:
   1. Die Funktion f(y) besitzt bei y₀ = 1 ein eindeutiges Minimum.
   2. Aufgrund der speziellen Gestalt des Integranden exp(-n·f(y) liefern die y-Werte aus der Umgebung von y₀ den Hauptbeitrag zum Integral.
   3. Daher wird die Funktion f(y) in ihre Taylor-Reihe entwickelt und durch die Approximation bis zur zweiten Ordnung erstetzt.
   4. Im dabei entstehenden Integral wird der Integrationsbereich auf alle reellen Zahlen erweitert. Es entsteht ein Gauß-Integral, das leicht ausgewertet werden kann.

Da im Endergebnis immer noch der Parameter n enthalten ist, hat man eine Funktion F(n) hergeleitet, die zur Approximation der Fakultät n! verwendet werden kann.

In Abbildung 1 wird versucht, die wichtigsten Rechenschritte und Zwischenergebnisse in einem Diagramm darzustellen; sie werden im Folgenden ausführlich erläutert und die Graphen der relevanten Funktionen werden gezeigt. Die blauen Pfeile in in Abbildung 1 stehen für Rechenschritte, bei denen eine Näherung gemacht wird, rote Pfeile symbolisieren Umformungen ohne Näherungen.

Abbildung 1: Übersicht über die Herleitung der Stirling-Approximation für die Fakultät mit Hilfe dr Laplace-Methode.

Die Berechnung der Stirling-Approximation

Die Gamma-Funktion

Die Gamma-Funktion wird als ein uneigentliches Integral nach Gleichung (1) in Abbildung 2 definiert. Man kann leicht zeigen, dass für x > 0 das Integral existiert. Die Eigenschaften der Gamma-Funktion sollen hier nicht ausführlich diskutiert werden; die Aufgaben unten enthalten einige Andeutungen.

Wichtig zur Berechnung der Stirling-Approximation der Fakultät sind folgende Eigenschaften der Gamma-Funktion:

1. Die Gamma-Funktion erfüllt die Rekursionsformel, die auch von der Fakultät erfüllt wird (siehe Gleichung (3) in Abbildung 2).
2. Setzt man in die Definition der Gamma-Funktion für x natürliche Zahlen ein, stimmen die uneigentlichen Integrale mit den Fakultäten überein (siehe Gleichung (4) in Abbildung 2).

Abbildung 2: Die Definition der Gamma-Funktion sowie ihre Eigenschaften, die benötigt werden, um die Stirling-Approximation für die Fakultät zu berechnen.

Aufgaben:

1. Zeigen Sie, dass die Gamma-Funktion die Rekursionsformel (3) in Abbildung 2 erfüllt.

2. Welche Integrale entstehen, wenn man die Gamma-Funktion nach Gleichung (1) explizit für natürliche Zahlen x = 1, 2, ... berechnet. Kann man diese Integrale explizit berechnen?

3. Welche Integrale entstehen, wenn man die Gamma-Funktion nach Gleichung (1) explizit für die Zahlen x = 1/2, 3/2, ... berechnet. Kann man diese Integrale explizit berechnen?

4. Skizzieren Sie die Integranden aus Gleichung (1) in Abbildung 2 für die Zahlen x aus der 2. und 3. Aufgabe.

5. Zeigen Sie dass man den Definitionsbereich der Gamma-Funktion in der Darstellung nach Gleichung (1) problemlos auf x > -1 erweitern kann.

Vorbereitungen zur Anwendung der Laplace-Methode

Der Ausgangspunkt für die Berechnung der Stirling-Approximation der Fakultät ist die Darstellung der Gamma-Funktion nach Gleichung (1) in Abbildung 2. Durch geeignete Substitutionen wird das Integral auf die Gestalt (9) in Abbildung 3 gebracht. Die Rechenschritte sind in Abbildung 3 dargestellt:

1. Gleichung (1) aus Abbildung 2 wird verwendet, um die Fakultät n! zu berechnen (siehe Gleichung (1) in Abbildung 3).
2. Im Faktor tⁿ wird die Basis t beseitigt (siehe Gleichung (1) in Abbildung 3).
3. Der Integrand wird als Exponentialfunktion geschrieben, wobei der Faktor exp(-n) vorgezogen wird (siehe Gleichungen (2) und (3) in Abbildung 3).
4. Durch die Substitution y = t/n nach Gleichung (4) wird das Integral in eine Integration über y verwandelt. Ziel der Substitution ist es, als Integranden eine Funktion in der folgenden Form zu erzeugen:

exp(-n·f(y))

mit einer geeigneten Funktion f(y), siehe Gleichung (9) in Abbildung 3. Die Rechenschritte (siehe Gleichung (5 - 8) in Abbildung 3) beinhalten die Rechengesetze für die Exponential- und Logarithmusfunktion.

Abbildung 3: Darstellung der Gamma-Funktion durch ein Integral (bezüglich y) über eine Funktion exp(-n·f(y) durch eine geeignete Substitution.

Die Laplace-Methode beruht darauf, das Integral in Gleichung (9) abzuschätzen, wobei die Funktion f(y) ein eindeutiges Minimum besitzen muss. Die folgenden Unterabschnitte diskutieren die Eigenschaften der Funktion f(y) und die Abschätzung des Integrals.

Die Eigenschaften der Funktion f(y) im Integranden

In Abbildung 3 wurde das uneigentliche Integral zur Berechnung der Gamma-Funktion mit Hilfe einer Funktion

f(y) = y - ln y

dargestellt. In Abbildung 4 werden ihre Eigenschaften untersucht. Für die Laplace-Methode ist entscheidend, dass f(y) ein eindeutiges Minimum besitzt.

Abbildung 4 zeigt den Nachweis, dass die Funktion f bei y₀ = 1 ein Minimum besitzt (siehe Gleichung (2 - 4)). Weiter wird die Taylor-Reihe bis zur zweiten Ordnung angegeben (bei Entwicklung um y₀ = 1; siehe Gleichung (5)).

Abbildung 4: Eigenschaften der Funktion f(y) = y - ln y. Für die Anwendung der Laplace-Methode ist wichtig, dass f(y) ein eindeutiges Minimum besitzt.

Abbildung 5 zeigt den Plot der Funktion f(y) = y - ln y (rot). Weiter eingetragen sind: die beiden Summanden, aus denen sich die Funktion f(y) zusammensetzt (die Ursprunsgerade mit Steigung 1 blau und der natürliche Logarithmus grün) sowie die Taylor-Reihe von f bis zur zweiten Ordnung (wobei um y₀ = 1 entwickelt wird, türkisfarben).

Abbildung 5: Plot der Funktion f(y), ihrer Summanden sowie der Taylor-Entwicklung bis zur zweiten Ordnung

Die näherungsweise Berechnung der Gamma-Funktion

Zuletzt wurde die Funktion f(y) in eine Taylor-Reihe entwickelt (um ihr Minimum y₀). Die eigentliche Laplace-Methode besteht jetzt darin, das Integral aus Gleichung (9) in Abbildung 3 näherungsweise zu berechnen. Sämtliche Rechenschritte sind in Abbildung 5 dargestellt.

1. Die Funktion f(y) wird durch ihre Taylor-Reihe (2. Ordnung) ersetzt (siehe Gleichung (1 - 4) in Abbildung 5).
2. Die Stelle y₀, bei der die Funktion f(y) ihr Minimum annimmt, befindet sich im Inneren des Integrationsbereiches. Der Integrand aus Gleichung (9) in Abbildung 3 besitzt daher bei y₀ ein Maximum.
3. Für sehr große n wird dieses Maximum sehr scharf und nur in einer kleinen Umgebung von y₀ entstehen die relevanten Beiträge zum Integral. Es sollte daher keinen großen Unterschied machen, wenn man den Integrationsbereich auf alle reellen Zahlen ausdehnt, siehe Gleichung (5). Durch eine weitere Substitution (siehe Gleichung (6) und (7) in Abbildung 5) man erhält ein Gaußintegral, das explizit berechnet werden kann, siehe Gleichung (7).

Abbildung 6: Umwandlung in ein Gauß-Integral und Berechnung des Gauß-Integrals.

Gleichung (9) zeigt jetzt das Endergebnis der Berechnung: die Stirling-Approximation für die Fakultät.

Aufgabe:

Schreiben Sie die rechte Seite der Stirling-Approximation (Gleichung (9) in Abbildung 6) als Exponentialfunktion (zur Basis e) und diskutieren Sie, wie die dadurch entstehende Funktion von n abhängt.

Veranschaulichung der relevanten Integranden

Oben wurde die Laplace-Methode zur Berechnung der Stirling-Approximation für die Fakultät gezeigt. Die Rechenschritte sollen nochmals erläutert werden; dazu werden die relevanten Integrale näher betrachtet und die Integranden dargestellt:

1. Der Integrand der Gamma-Funktion: siehe Gleichung (1) in Abbildung 7 und die Graphen in Abbildung 8.
2. Der Integranden aus Gleichung (9) in Abbildung 3: Er entsteht, wenn im Integral der Gamma-Funktion eine Substitution durchgeführt wird: Die Substitution wird dabei so gewählt, dass ein Integrand entsteht, so dass nahezu alle Beiträge zum Integral aus einem kleinen Bereich um das Minimum der Funktion f(x) = x - ln x stammen. Siehe Gleichung (2) in Abbildung 7 und die Graphen in Abbildung 9.
3. Der Integrand aus Gleichung (5) in Abbildung 6. Jetzt wurde die Funktion f(x) = x - ln x in eine Taylor-Reihe entwickelt: siehe Gleichung (3) in Abbildung 7 und die Graphen in Abbildung 10.
4. Der Integrand aus Gleichung (5) in Abbildung 6. Dazu wurde ein Gauß-Integral erzeugt, indem der Integrand verschoben (er ist jetzt symmetrisch zur y-Achse) und skaliert wurde (in der von x abhängigen Exponentialfunktion tritt kein n mehr auf). siehe Gleichung (4) in Abbildung 7 und die Graphen in Abbildung 11.

Abbildung 7: Die relevanten Integranden aus den Integralen der Gamma-Funktion zur Berechnung von n! beziehungsweise der Näherungen gemäß der Laplace-Methode für n = 1, 2,..., 8.

Es wurde mehrmals betont, dass die Stirling-Approximation verwendet wird, um die Fakultät n! für große n zu berechnen. In den folgenden Abbildungen werden aber sehr kleine n eingesetzt (n = 1, 2, ..., 8). Aber selbst hier lässt sich nachvollziehen, wie die Laplace-Methode arbeitet und wie ein Integral entsteht, so dass die zu berechnende Fläche in einem kleinen Intervall konzentriert ist.

Der Integrand der Gamma-Funktion

Zur Berechnung der Fakultäten

n! für n = 1, 2,..., 8

mit Hilfe der Gamma-Funktion Γ(n+1) Gleichung (1) aus Abbildung 2 müssen die Integranden

g_n(x) = xⁿ·e^-x, für n = 1, 2,..., 8

betrachtet werden (jetzt mit x als Variable, siehe auch Gleichung (1) in Abbildung 7). Diese Integranden sind in Abbildung 8 dargestellt. Die Funktion xⁿ ist streng monoton ansteigend und geht gegen unendlich für x → ∞. Dagegen ist die Funktion e^-x streng monoton abnehmend und geht gegen null für x → ∞. Als Produkt entsteht eine Funktion mit einem eindeutigen Maximum (Ausnahme ist hier n = 1). Die Fläche unter der Kurve entspricht der Fakultät n!. Mit "Fläche" ist hier das uneigentliche Integral gemeint, wenn die Integrationsvariable die Werte von 0 bis ∞ durchläuft.

Abbildung 8: Integranden der Gamma-Funktion zur Berechnung von n! für n = 1, 2,..., 8.

Man erkennt, dass das Maximum von g_n(x) für zunehmende n sehr breit wird wird und das Intervall, das den Großteil zum Integral beiträgt, mit n anwächst. Schwerer zu erkennen ist (aufgrund der unterschiedlichen Skalierungen), dass die Funktionswerte beim Maximum mit n stark ansteigen.

Aus den Abbildungen lässt sich natürlich nicht ersehen, wie die Fakultäten n! von n abhängen. Aber der Vergleich der Funktionen aus Abbildung 8 mit den Integranden nach der Substitution beziehungsweise Einsetzen der Taylor-Reihe wird die Laplace-Methode besser veranschaulichen.

Der Integrand der Gamma-Funktion nach der Substitution

In Abbildung 3 wurde gezeigt, wie das Integral der Gamma-Funktion durch eine Substitution umgeformt wurde (siehe Gleichung (9) in Abbildung 3); jetzt lautet der Integrand zur Berechnung von Γ(n+1) = n! (wieder als Funktion von x; siehe auch Gleichung (2) in Abbildung 7):

f_n(x) = nⁿ⁺¹·exp(-n(x - ln x), für n = 1, 2,..., 8.

In Abbildung 9 werden diese Integranden dargestellt.

Abbildung 9: Integranden der Gamma-Funktion zur Berechnung von n! für n = 1, 2,..., 8 nach der Substitution.

Die Abbildungen sind wieder schwer miteinander zu vergleichen, da die Skalierungen an die Funktionswerte angepasst wurde. Aber im Vergleich zu Abbildung 8 erkennt man:

1. Die Funktionen besitzen immer noch ein eindeutiges Maximum. Jetzt aber unabhängig von n immer bei x = 1.
2. Die Funktionswerte beim Maximum wachsen schnell mit n an.
3. Aber die Funktionswerte, die den größten Beitrag zum Integral liefern, sind sehr nahe um das Maximum bei x = 1 konzentriert. Qualitativ besitzen die Integranden eine nahezu identische Breite (wie immer man diese "Breite" definiert).

Es soll aber nochmals betont werden: Bisher wurde keine Näherung verwendet. Durch die Substitution wird das Integral lediglich transformiert, der Wert des Integrals liefert immer noch den exakten Wert für die Fakultät.

Näherung: Der Integrand nach dem Einsetzen der Taylor-Entwicklung

Die Integranden aus Abbildung 9 werden jetzt approximiert. Dazu wird die Funktion x - ln x um ihr Minimum (bei x = 1) in die Taylor-Reihe bis zur zweiten Ordnung entwickelt (zur Veranschaulichung: siehe Abbildung 5) und in das Argument der Exponentialfunktion eingesetzt.

Das heißt die Funktion

f(x) = x - ln x

wird bei x = 1 in eine Taylor-Reihe entwickelt, wodurch eine quadratische Funktion entsteht:

f(x) = x - ln x ≈ 1 + (x - 1)²/2.

Anstelle der Integranden f_n(x) erhält man jetzt (siehe Gleichung (3) in Abbildung 7):

h_n(x) = nⁿ⁺¹·exp(-n)·exp(-n·(x - 1)²/2), für n = 1, 2,..., 8.

Die Funktionen h_n(x) werden in Abbildung 10 dargestellt.

Abbildung 10: Integranden nach Einsetzen der Taylor-Reihe zur Berechnung von n! für n = 1, 2,..., 8.

Vergleicht man diese Graphen mit denen aus Abbildung 9, so erkennt man:

1. Der Definitionsbereich ist jetzt nicht mehr auf die positiven reellen Zahlen beschränkt, sondern kann auf die gesamten reellen Zahlen ausgedehnt werden.
2. Die Integranden in Abbildung 9 und Abbildung 10 besitzen jeweils ein eindeutiges Maximum. Da die Taylor-Entwicklung um das Minimum von f(x) = x - ln x zur Näherung verwendet wurde, ist die Lage des Maximums jeweils unverändert.
3. Aber da die Taylor-Entwicklung bis zur zweiten Ordnung eine quadratische Funktion ist, entstehen in Abbildung 10 Funktionen, die zu x = 1 achsensymmetrisch sind.
4. Da die einzelnen Diagramme wieder unterschiedlich skaliert sind, ist es schwer die einzelnen Graphen miteinander zu vergleichen. Bei genauerer Betrachtung erkennt man aber, dass mit zunehmendem n die Graphen schlanker werden und somit der Hauptbeitrag zum Integral über die Funktion aus einem kleinen Bereich um das Maximum stammt.
5. Wie gut die Approximation der Integranden aus Abbildung 9 nach dem Einsetzen der Taylor-Entwicklung ist, lässt sich hier nur schwer ablesen; einen besseren Vergleich werden die Abbildungen weiter unten liefern, in denen dann die Graphen aus Abbildung 9 und 10 gemeinsam dargestellt werden.

Die Funktionen h_n(x) beinhalten eine Exponentialfunktion, in der sowohl der Parameter n als auch die Variable x vorkommen. Die Berechnung der Integrale über die Funktionen h_n(x) wird übersichtlicher, wenn man stattdessen eine Produktdarstellung finden kann, in der ein Faktor nur von n und der andere nur von x abhängt. Dazu wird eine weitere Substitution durchgeführt (siehe Gleichung (6) in Abbildung 6), wodurch sich der von n abhängige Faktor komplett vor das Integral ziehen lässt und ein von n unabhängiges Gauß-Integral entsteht. Durch die Substitution entstehen Funktionen k_n(x), die in Gleichung (7) in Abbildung 6 und Gleichung (4) in Abbildung 7 gezeigt sind:

k_n(x) = n^n+1/2·exp(-n)·exp(-x²/2), für n = 1, 2,..., 8.

Gauß-Integrale werden hier nicht besprochen; ihre Berechnung wird in Eigenschaften von Zufallsvariablen: Die Varianz und die Standardabweichung ausführlich gezeigt.

Abbildung 11 zeigt die Integranden nachdem die Funktionen h_n(x) in Gaußsche Glockenkurven verwandelt wurden, die zu x = 0 achsensymmetrisch sind.

Abbildung 11: Integranden nach Einsetzen der Taylor-Reihe zur Berechnung von n! für n = 1, 2,..., 8 und einer weiteren Substitution, die Gaußsche Glockenkurven erzeugt, die zu x = 0 achsensymmetrisch sind.

Aufgabe: Die Graphen in den Abbildungen 9 und 10 besitzen offensichtlich jeweils zwei Wendepunkte. Berechnen Sie die Koordinaten der Wendepunkte.

Kann man die "Breite" der Funktionen durch den Abstand der Wendepunkte in x-Richtung definieren? Wie hängt diese Breite für die Funktionen h_n(x) und k_n(x) von n ab? Gegen welchen Wert geht diese Breite für n → ∞?

Vergleich der Integranden vor und nach der Taylor-Entwicklung

Der entscheidende Schritt in der Berechnung der Stirling-Approximation für die Fakultät ist somit die Taylor-Entwicklung der Funktion

f(x) = x - ln x

bis zur zweiten Ordnung, denn er ermöglicht es, das Integral der Gamma-Funktion durch ein Gauß-Integral zu approximieren. Die Integranden vor und nach der Taylor-Entwicklung sind in Abbildung 9 und 10 dargestellt, also die Graphen der Funktionen f_n(x) und h_n(x).

Abbildung 12 zeigt jetzt für die Werte n = 1, 2,..., 8 die Funktionen f_n(x) (blau) und h_n(x) (dunkelrot) in jeweils einem Diagramm.

Abbildung 12: Die Integranden f_n(x) (blau) und h_n(x) (rot) für n = 1, 2,..., 8 vor und nach der Taylor-Entwicklung von f(x) = x - ln x.

Man erkennt, dass für kleine Werte von n noch deutliche Unterschiede betsehen. Aber schon bei n = 8 sind in der gezeigten Skalierung die Graphen nur bei kleinen Funktionswerten zu unterscheiden. Insbesondere erkennt man, dass die "Breite" der Funktionen (also der Abstand der Wendepunkte in x-Richtung) nahezu identisch sein muss.

In Abbildung 13 und 14 werden dann die entsprechenden Graphen für n = 20 beziehungsweise n = 100 gezeigt. Zusätzlich ist die Gaußsche Glockenkurve gezeigt (türkisfarben), die bei der letzten Substitution entsteht (Übergang von h_n(x) zu k_n(x)).

Abbildung 13: Die Integranden aus Abbildung 12 für n = 20. Zusätzlich ist der Integrand nach der letzten Substitution gezeigt (türkisfarben).

Abbildung 14: Die Integranden aus Abbildung 12 für n = 100. Zusätzlich ist der Integrand nach der letzten Substitution gezeigt (türkisfarben).