Verwendung von Fourier-Koeffizienten zum Testen auf Gleichverteilung

Jan Behrens
1. Juni 2008
(Ergänzung/Überarbeitung vom 12. Mai 2010)

 

Es wird unentgeltlich jeder Person das Recht eingeräumt, Kopien dieses Dokumentes und des dazugehörigen Haskell-Programmes zu benutzen, zu kopieren, zu modifizieren, mit anderen Werken zu kombinieren, zu veröffentlichen, zu verbreiten, unterzulizensieren, zu verkaufen und/oder weiteren Personen diese Rechte einzuräumen, soweit in allen Kopien oder wesentlichen Teilen davon der Autor genannt wird und dieser Hinweistext enthalten bleibt. Fehlerfreiheit und Funktionsfähigkeit werden nicht zugesichert; Benutzung erfolgt auf eigenes Risiko.

Dieses Dokument ist auch als PDF-File (10 Seiten, DIN-A4) verfügbar, außerdem kann die LaTeX-Quelldatei heruntergeladen werden.


 

Ziel soll es sein, zu einer Stichprobe von Zufallszahlen im Bereich von 0 bis 1 den Wert g(h) ∈ [0,1] zu berechnen, der als Maß für das Abweichen dieser Stichprobe von einer Gleichverteilung dient. Kleine Werte von g(h) bedeuten eine hohe Abweichung von der Gleichverteilung, während Werte gegen 1 für eine gute Übereinstimmung der Stichprobe mit einer Gleichverteilung stehen. Die Größe g(h) soll außerdem bei größeren Stichproben (M > 30) näherungsweise das Signifikanzniveau wiederspiegeln, d.h. die Wahrscheinlichkeit, dass g(h) bei zu Grunde liegender Gleichverteilung einen Wert x unterschreitet, soll näherungsweise x betragen.

Hierzu wird die Summe h der Quadrate von Fourier-Koeffizienten (unter Auslassung des Koeffizienten für den Gleichanteil) gebildet, da diese bei einer konstanten Funktion (entsprechend einer Gleichverteilung) gleich 0 ist. Konvergenz wird dadurch gewährleistet, dass mit zunehmendem Index n die Quadrate der Fourier-Koeffizienten schwächer gewichtet werden, nämlich mit 1/n². Zur Berechnung der Größe g(h) aus der Summe h wird diese Summe für den Fall einer zugrundeliegenden Gleichverteilung (hgleich) zunächst mit Hilfe des Zentralen Grenzwertsatzes von Lindeberg-Levy durch eine unendliche Summe von gewichteten, quadrierten normalverteilten Zufallsgrößen angenähert. Die der unendlichen Summe dieser Zufallsgrößen entsprechenden unendlichen Faltungen der Wahrscheinlichkeitsdichten werden anschließend berechnet, um die Wahrscheinlichkeitsdichte f(x) von hgleich zu ermitteln. Aus dieser Wahrscheinlichkeitsdichte wird anschließend durch Integration die kumulative Verteilungsfunktion F(x) und daraus dann eine Formel für die gesuchte Größe g(h) gebildet.

Es wird weiterhin eine Möglichkeit erwähnt, wie durch eine einfache Transformation der Stichprobe anstelle des Testens auf Gleichverteilung auf jede andere Wahrscheinlichkeitsverteilung mit gegebenen Parametern geprüft werden kann.

Da die Signifikanz bei Erhöhung des Stichproblenumfanges von einer nicht exakt gleichverteilten Zufallsgröße beliebig steigt und g(h) entsprechend gegen 0 strebt, werden zum Abschluss noch weitere Größen k und k* untersucht, welche ebenfalls Maße für die Ungleichverteilung sind, sich jedoch bei einer Erhöhung des Stichproblemumfanges einem Grenzwert annähern, der die Abweichung zwischen der zugrundeliegenden Verteilung und einer Gleichverteilung wiederspiegelt.

 

Gegeben sei eine Folge a von M reellen Zahlen, die im Intervall [0;1] liegen. Durch Berechnung von

h := 1/M * SUM n=1..inf [1/n^2 * ((SUM m=1..M cos(2pi*n*a[m]))^2 + (SUM m=1..M sin(2pi*n*a[m]))^2]

erhalten wir eine Zahl, die umso größer wird, desto ungleich-verteilter die Zahlen der Folge a im Intervall [0;1] sind, da bei einer Gleichverteilung sich die positiven und negative Werte des Sinus und Cosinus in den (inneren) Summen aufheben.

Sind die Zahlen der Folge zufällig und genügen einer Gleichverteilung, dann weisen die Zufallsgrößen

U[n,m] := cos(2pi*n*a[m])

folgenden Mittelwert μ und folgende Varianz σ² auf:

mu = 0
sigma = 1/2

Wir bilden nun die Zufallsgröße:

V[n] := (SUM m=1..inf U[n,m]) / sqrt(M) = SUM m=1..M cos(2pi*n*a[m]) / sqrt(M)

Vn nähert sich für M → ∞ gemäß dem Grenzwertsatz von Lindeberg-Levy einer Normalverteilung mit dem Mittelwert μ = 0 und der Varianz σ² = ½ an. Für M > 30 erhalten wir mit der Normalverteilung bereits eine hinreichend gute Näherung. Seien Xn und Yn (0,1)-normalverteilte Zufallsgrößen, dann läßt sich Vn (näherungsweise) auch darstellen als:

V[n] := mu + sigma*X[n] = sqrt(1/2) * X[n]

Die Zufallsgrößen

W[n] := 1/n^2 * ((SUM m=1..M cos(2pi*n*a[m]) / sqrt(M))^2 + (SUM m=1..M sin(2pi*n*a[m]) / sqrt(M))^2)
= 1/M * SUM n=1..inf [1/n^2 * ((SUM m=1..M cos(2pi*n*a[m]))^2 + (SUM m=1..M sin(2pi*n*a[m]))^2]

ergeben sich entsprechend zu:

W[n] = 1/n^2 * ((sqrt(1/2) * X[n])^2 + (sqrt(1/2) * Y[n])^2) = 1/n^2 * (X[n]^2 + Y[n]^2) / 2

h lässt sich somit bei gleichverteilt-zufälligen am als unendliche Summe von gewichteten, quadrierten (0,1)-normalverteilten Zufallsgrößen Xn und Yn darstellen:

h[gleich] = SUM n=1..inf W[n] = SUM n=1..inf 1/n^2 * (X[n]^2 + Y[n]^2) / 2

Die Zufallsgröße Xn² + Yn² genügt einer χ²-Verteilung mit 2 Freiheitsgraden, welche die Wahrscheinlichkeitsdichte ½ · exp(−½ · x) aufweist. Multipliziert man eine Zufallsgröße mit einem Faktor, so bedeutet dies für die Wahrscheinlichkeitsdichte eine entsprechende Streckung in x-Richtung sowie eine Stauchung in y-Richtung. Demnach entsprechen die einzelnen Summanden

1/n^2 * (X[n]^2 + Y[n]^2) / 2

von hgleich den folgenden Wahrscheinlichkeitsdichtefunktionen:

n^2 * exp(-n^2 * x)

Ein Addieren von Zufallsgrößen entspricht der Faltung (∗) ihrer Wahrscheinlichkeitsdichten. Bezeichnen wir mit f(x) die Wahrscheinlichkeitsdichte der Zufallsgröße hgleich an der Stelle x, so ist:

f(x) = (exp(-x)) (*) (2^2 * exp(-2^2 * x) (*) (3^2 * exp(-3^2 * x) (*) (4^2 * exp(-4^2 * x) (*) ...

mit
r(x) (*) s(x) = INTEGRAL Theta=0..x r(x-Theta) * s(Theta) dTheta

Wir berechnen zunächst die Faltung von zwei unterschiedlichen Exponentialfunktionen (αβ):

(A * exp(-alpha * x)) (*) (B * exp(-beta * x)) = [...]
[...]
[...] = (A * B) / (alpha - beta) * (exp(-beta * x) - exp(-alpha * x))

Wir stellen fest, dass die Faltung zweier Exponentialfunktionen mit unterschiedlich skaliertem Exponenten eine Linearkombination dieser Exponentialfunktionen ergibt.

Faltet man eine Summe von Exponentialfunktionen mit einer weiteren Exponentialfunktion A · exp(−αx), so ist der Koeffizient B eines Summanden B · exp(−βx) im Ergebnis mit A / (αβ) zu multiplizieren. Aus diesem Grunde lässt sich f(x) für x > 0 darstellen als:

f(x) = SUM n=1..inf (PRODUCT m element N\{n} (m^2)/(m^2-n^2)) * (n^2 * exp(-n^2 * x))

Zur Berechnung des unendlichen Produktes p(n) spalten wir dieses zunächst in drei getrennte Produkte p1(n), p2(n) und p3(n) auf:

p(n) = PRODUCT m element N\{n} (m^2)/(m^2-n^2) = PRODUCT m=1..n-1 [...] * PRODUCT m=n+1..2n [...] * PRODUCT m>2n [...]

Umformung von p1(n) ergibt:

p1(n) = [...]
[...]
[...] = (-1)^(n-1) * PRODUCT m=1..n-1 m/(m+n)

Umformung von p2(n) ergibt:

p2(n) = [...] = PRODUCT m=1..n ((m+n)^2)/(m*(m+2n))

Und eine Umformung des unendlichen Produktes p3(n) ergibt:

p3(n) = [...]
[...]
[...] = PRODUCT m=1..n (m+2n)/(m+n)

Das Gesamtprodukt p(n) ist somit:

p(n) = [...]
[...]
= (-1)^(n-1) * 2

Setzen wir dieses Ergebnis in die Gleichung für f(x) ein, so erhalten wir (x > 0 vorausgesetzt):

f(x) = SUM n=1..inf p(n) * (n^2 * exp(-n^2 * x))
[...] = 2 * SUM (-1)^(n-1) * n^2 * exp(-n^2 * x)

Aus der Dichteverteilung f lässt sich die kumulative Verteilungsfunktion F durch Integration berechnen:

F(x) = [...] = 1 - 2 * SUM n=1..inf (-1)^(n-1) * exp(-n^2 * x) für x > 0, sonst 0

Mit der gegebenen Verteilungsfunktion F lässt sich nun leicht die Wahrscheinlichkeit g(x) berechnen mit der die Größe h für den Fall, dass die Zahlen der Folge am die Stichprobe einer gleichverteilten Zufallsgröße sind, einen gegebenen Wert x überschreitet:

g(x) = P(h[gleich] > x) = 1 - F(x) = 2 * SUM n=1..inf (-1)^(n-1) * exp(-n^2 * x) für x > 0, sonst 1

Ist g(h) < 5%, dann weicht die Stichprobe signifikant von einer Gleichverteilung ab, d.h. die Wahrscheinlichkeit beträgt nur 5%, dass im Falle des Zugrundeliegens einer Gleichverteilung h zufällig so groß bzw. g(h) zufällig so klein ist.

Das beschriebene Verfahren kann nicht nur zur Prüfung einer Stichprobe auf Gleichverteilung, sondern ebenso zur Prüfung einer Stichprobe a′ auf jede beliebige Verteilung (mit gegebenen Parametern) angewendet werden. Hierzu ist die Stichprobe vorher mit Hilfe der entsprechenden kumulativen Verteilungsfunktion Ψ auf die geprüft werden soll zu transformieren:

a[m] = Psi(a'[m])

Möchte man auf eine Normalverteilung mit dem Mittelwert μ und der Standardabweichung σ prüfen, dann ist:

a[m] = Psi(a'[m]) = Phi[mu,sigma](a'[m]) = 1/2 * (1 + erf((a'[m] - mu) / (sigma * sqrt(2))))

h bzw. g(h) eignen sich, um festzustellen wie signifikant eine Stichprobe von einer Gleichverteilung abweicht. Da, falls die Stichprobe keiner Gleichverteilung folgt, die Signifikanz bei Erhöhung des Stichprobenumfangs beliebig steigt, definieren wir noch ein Maß k ∈ [0,1] für den Grad der Abweichung von einer Gleichverteilung, welches sich bei Erhöhung des Stichprobenumfanges einem Grenzwert annähert:

k := sqrt((6 / pi^2 / M^2) * M * h)

Insbesondere wirkt sich ein Duplizieren der Stichprobe nicht auf die neu definierte Größe k aus, d.h. für am+M0 = am und NN gilt k|M = N·M0 = k|M = M0. Falls eine Stichprobe mit M optimal gleichverteilten Werten vorliegt (am = m/M + d mit d ∈ [0,m/M]), ist die Größe k = 1/M. Dies lässt sich wie folgt beweisen:

k|a[m]=m/M+d = sqrt((6 / pi^2 / M^2) * SUM n=1..inf (1/n^2 * k[m]))
k[m] := (SUM n=1..M cos(2*pi*n*(m/M+d)))^2 + (SUM n=1..M sin(2*pi*n*(m/M+d)))^2

Man formt den Term km mittels Ausmultiplizieren und Additionstheoremen um:

k[m] = SUM m=1..M SUM m'=1..M [ cos(2*pi*n*(m/M+d)) * cos(2*pi*n*(m'/M+d)) ] + SUM m=1..M SUM m'=1..M [ sin(2*pi*n*(m/M+d)) * sin(2*pi*n*(m'/M+d)) ]
= SUM m=1..M SUM m'=1..M [ cos(2*pi*n*(m/M+d)) * cos(2*pi*n*(m'/M+d)) + sin(2 *pi*n*(m/M+d)) * sin(2*pi*n*(m'/M+d)) ]
= SUM m=1..M SUM m'=1..M cos(2*pi*n*(m/M+d) - 2*pi*n*(m'/M+d))

Dann lässt sich die Variable d entfernen:

= SUM m=1..M SUM m'=1..M cos(2*pi*n*m/M - 2*pi*n*m'/M)

Anschließend wird der Term zurück transformiert:

k[m] = SUM m=1..M SUM m'=1..M [ cos(2*pi*n*m/M) * cos(2*pi*n*m'/M) + sin(2 *pi*n*m/M) * sin(2*pi*n*m'/M) ]
= SUM m=1..M SUM m'=1..M [ cos(2*pi*n*m/M) * cos(2*pi*n*m'/M) ] + SUM m=1..M SUM m'=1..M [ sin(2*pi*n*m/M) * sin(2*pi*n*m'/M) ]
= (SUM n=1..M cos(2*pi*n*m/M))^2 + (SUM n=1..M sin(2*pi*n*m/M))^2

Falls n ein Vielfaches von M ist, ist km = M², ansonsten sind beide Summen 0 und damit auch km = 0. Beschreiben wir die Vielfachen von M mit i·M, dann ergibt sich für k|am = m/M + d:

k|a[m]=m/M+d = sqrt((6 / pi^2 / M^2) * SUM i=1..inf (1/(i*M)^2 * M^2))
= sqrt((6 / pi^2 / M^2) * SUM i=1..inf 1/i) = 1/M

Dieses entspricht genau der oben aufgestellten Behauptung.

Wie gezeigt wurde, nimmt die Größe k bei perfekter Gleichverteilung einer Stichprobe den Wert 1/M an. Kleinere Werte sind bei einem endlichen Stichprobenumfang M nicht möglich. Wir definieren eine weitere Größe k* deren Werte ebenfalls im Intervall [0,1] liegen:

k* := sqrt((M * h * 6 / pi^2 - 1) / (M^2 - 1))

Die Größe k* nimmt im Falle einer optimalen Gleichverteilung der Stichprobe im Gegensatz zur Größe k nicht den Wert 1/M sondern den Wert 0 an. Es gilt also: k*|am = m/M + d = 0. Der Beweis hierfür ist im Wesentlichen identisch mit dem Beweis, dass k|am = m / M + d = 1/M ist. Für große M ist kk*.


Folgendes Haskell-Programm (hier herunterzuladen) kann die Größen g(h) und k* berechnen:

inhomoSum :: RealFloat a => [a] -> a
inhomoSum xs
  =  sum [ ((s1 n)**2 + (s2 n)**2) / n**2
         | i <- [1..huge], n <- [fromIntegral i] ]
  where
    s1 n  =  sum [ cos (2 * pi * n * x) | x <- xs ]
    s2 n  =  sum [ sin (2 * pi * n * x) | x <- xs ]
    huge  =  1000

inhomoSignificance :: RealFloat a => [a] -> a
inhomoSignificance xs
  | xs == []   =  error "inhomoSignificance undefined for empty list"
  | otherwise  =  prob (inhomoSum xs / count)
  where
    count   =  fromIntegral (length xs)
    prob x  =  limit (
                 sum [ (if i < huge then 2 else 1) *
                       (-1)**(n-1) * exp (-n**2 * x)
                     | i <- [1..huge], n <- [fromIntegral i] ]
               )
    limit   = (min 1 . max 0)
    huge    =  1000

inhomogeneity :: RealFloat a => [a] -> a
inhomogeneity xs
  | xs == []        =  error "inhomogenity undefined for empty list"
  | length xs == 1  =  error "inhomogenity undefined for single value"
  | otherwise       =  sqrt ( limit (
                         (inhomoSum xs * 6 / pi**2 - 1) /
                         (count**2 - 1)
                       ) )
  where
    count  =  fromIntegral (length xs)
    limit  =  max 0

Hinweise zum Programm: fromIntegral dient nur der Typ-Konvertierung von ganzzahligen Zahlen in Gleitkommazahlen und hat nichts mit der Berechnung von mathematischen Integralen zu tun. Die Variablen i und n stellen beide die gleiche Zahl dar, jedoch unterschiedlichen Typs: i ist vom Typ Integer, während n einen Typ der Klasse RealFloat aufweist. Die Fallunterscheidung  if i < huge then 2 else 1  dient der Vermeidung von Berechnungsfehlern für x gleich oder nahe 0, indem das letzte Glied der (eigentlich unendlichen) Summe nur halb angerechnet wird. Die Funktionen namens limit sind Hilfsfunktionen um numerische Ungenauigkeiten am Rande des Definitionsbereiches aufzulösen.

inhomoSum berechnet eine Zwischengröße, welche von den beiden anderen Funktionen inhomoSignificance und inhomogeneity verwendet wird. Die Funktion inhomoSignificance berechnet die Größe g(h) für eine gegebene Zahlenfolge, d.h. wenn das Ergebnis kleiner als 0.05 ist, liegt eine signifikante Abweichung von der Gleichverteilung vor. Die Funktion inhomogeneity berechnet die Größe k*, wobei 1 für maximale Ungleichverteilung (alle am sind identisch) und 0 für optimale Gleichverteilung (am = m/M + d) steht.


Zurück zur Hauptseite

Impressum