Forum:  Programmieren
Thema: Python: Cumulative Distribution Function implementieren
Themen-Übersicht
Lucky_7
Aktiv
Dabei seit: 21.01.2018
Mitteilungen: 143
Aus:
Themenstart: 2019-04-26 11:38

Hey,

ich versuche gerade einen so genannten Kolmogorov-Smirnov Test zu implementieren.

Hierbei habe ich eine geschätze Wahrscheinlichkeitsdichte-Funktion (PDF) gegeben und möchte die Cumulative-Distribution-Function (CDF) hieraus bilden. Das sollte eigentlich nur das Integral der PDF sein.
Ich habe gedacht, das würde ich durch diese Zeile:
python
cdf[idx] = np.sum(pdf[:idx])
im unten stehenden Code erreichen.

Irgendetwas mache ich aber falsch, weil meine CDF, die ich aus der geschätzten PDF berechne, Werte deutlich größer 1 annimmt:


Ich erkenne wirklich nicht, wo der Fehler liegt.
Dies ist der Code, den ich verwende:
python
cdf = np.zeros(len(pdf))
hypoCdf = np.zeros(len(pdf))
d_u = np.zeros(len(pdf)) #upper
d_l = np.zeros(len(pdf)) #lower
alpha = 0.05
for idx in range(0, len(pdf)):
    cdf[idx] = np.sum(pdf[:idx])
    hypoCdf[idx] = norm(loc=mean(pdf), scale=np.std(pdf)).cdf(idx) #cdf of gaussian distribution that we assume as our null hypothesis 
for idx in range(0, len(cdf)):
    d_u[idx] = abs(cdf[idx] - hypoCdf[idx])
    if idx == 0:
        d_l[idx] = 0  # is this correct ? 
    else:
        d_l[idx] = abs(cdf[idx-1] - hypoCdf[idx])
d_crit = np.sqrt(-0.5*np.log(alpha/2))/np.sqrt(len(pdf))
d_u_max = np.max(d_u)
d_l_max = np.max(d_l)
if d_u_max > d_crit or d_l_max > d_crit:
    print("result: null hypothesis rejected")
else:
    print("result: failed to reject null hypothesis")
plt.figure()
plt.plot(x_grid, cdf, linewidth = 1, alpha = 1, color = 'r')
plt.plot(x_grid, hypoCdf, linewidth = 1, alpha = 1, color = 'k')
plt.show()

..wenn euch das komplette Skript, oder die Daten interessieren, die ich zum Erzeugen der pdf verwendet habe, findet ihr einen Link in dieser Frage hier


Lucky_7
Aktiv
Dabei seit: 21.01.2018
Mitteilungen: 143
Aus:
Beitrag No.1, vom Themenstarter, eingetragen 2019-04-26 17:20

ok, ich denke es ist notwendig die gebildete Summe noch durch die Gesamtsumme zu teilen:
python
for idx in range(0, len(pdf)):
    if idx == 0:
        cdf[idx] = pdf[idx]/np.sum(pdf)
    else:
        cdf[idx] = np.sum(pdf[:idx])/np.sum(pdf)
    hypoCdf[idx] = norm(loc=mean(pdf), scale=np.std(pdf)).cdf(x_grid[idx])

Davon abgesehen, habe ich auch die hypothetische CDF falsch ausgewertet.
Für die CDF der geschätzten PDF dürfte es so aber stimmen. ..


schnitzel
Aktiv
Dabei seit: 26.02.2009
Mitteilungen: 140
Aus:
Beitrag No.2, eingetragen 2019-04-26 19:33

Hi,

Schleifen über Indices ist in python oft vermeidbar. Wenn man Schleifen benutzt, dann am besten direkt über die Elemente.
Hier sollte das aber auch mit der Funktion np.cumsum funktionieren.

Gruß


Lucky_7
Aktiv
Dabei seit: 21.01.2018
Mitteilungen: 143
Aus:
Beitrag No.3, vom Themenstarter, eingetragen 2019-04-27 12:36

Ok. Um d_l zu berechnen, benötige ich aber weiterhin eine Schleife - sehe ich das richtig ?
python
alpha = 0.05 # 5% significance level
cdf = np.cumsum(pdf)/np.sum(pdf)
hypoCdf = norm(loc=mean(pdf), scale=np.std(pdf)).cdf(x_grid) #cdf of gaussian distribution that we assume as our null hypothesis 
d_u = abs(cdf - hypoCdf) # upper limit
d_l = np.zeros(len(cdf))
for idx in range(0, len(cdf)):
    if idx == 0:
        d_l[idx] = 0  
    else:
        d_l[idx] = abs(cdf[idx-1] - hypoCdf[idx]) #lower limit
d_crit = np.sqrt(-0.5*np.log(alpha/2))/np.sqrt(len(pdf))
d_u_max = np.max(d_u)
d_l_max = np.max(d_l)
if d_u_max > d_crit or d_l_max > d_crit:
    print("result: null hypothesis rejected")
else:
    print("result: failed to reject null hypothesis")

... kannst du mir erklären, warum der KS-Test überhaupt d_l berechnet? Sollte es nicht eigentlich genügen, sich d_u anzuschauen? d_u vergleicht doch die CDF der Hypothese und der Sample-Distribution. Wozu gibt es dann noch d_l ?


b_p
Senior
Dabei seit: 11.09.2012
Mitteilungen: 541
Aus: Köln
Beitrag No.4, eingetragen 2019-07-01 00:35

Kurze Anmerkung am Rande: pdf scheint eine Liste bzw. ein Array mit den Beobachtungen zu sein. Die hypothetische Verteilung soll eine Normalverteilung sein, wobei deren Parameter offenbar nicht bekannt sind. Letztere werden folglich einfach aus den Beobachtungen geschätzt und dann für die hypothetische Verteilung verwendet.

Betrachten wir mal ganz allgemein die zwei grundsätzlich möglichen Fälle bei einem Anpassungstest:

Im ersten Fall sind sowohl die Verteilungsfamilie als auch ihre Parameter bekannt, d. h. die Nullhypothese ist vollständig spezifiziert und somit einfach. Dieser Fall ist in der Regel gut untersucht und die Verteilung der Teststatistik (zumindest asymptotisch) bekannt. Folglich entspricht dies auch dem in Lehre und Standardliteratur am häufigsten dargestellten Fall.

Im zweiten Fall, welcher eher der Praxis entspricht, liegt zwar das Wissen um die Verteilungsfamilie, nicht aber um deren Parameter vor. Die Nullhypothese ist folglich unvollständig spezifiziert und somit zusammengesetzt. Die Verteilung der Teststatistik weicht in der Regel erheblich von der aus dem vollständig spezifizierten Fall ab, wodurch das Festhalten an der alten Prüfverteilung in einem wesentlich geringeren, tatsächlichen Fehler 1. Art resultieren würde (vgl. S. 73f in Büning, H. und G. Trenkler (1994). Nichtparametrische statistische Methoden (2., erweiterte und völlig überarbeitete Auflage). Berlin, de Gruyter). Die korrekte Prüfverteilung hängt vom Stichprobenumfang, der hypothetischen Verteilung, den wahren Parameterwerten und der verwendeten Schätzmethode ab (vgl. S. 102 in D'Agostino, R. B. und M. A. Stephens (1986). Goodness-of-Fit Techniques. New York, Marcel Dekker). Handelt es sich um Lage- oder Skalenparameter, so werden bei einer Schätzung durch die "richtige" Methode zumindest die beiden zuletzt genannten Abhängigkeit eliminiert. Oftmals ist die Prüfverteilung jedoch nicht einmal bekannt, so dass Quantile und p-Werte mühsam über MC-Simulationen berechnet werden müssen.


Aus obigem wird klar, dass hier leider der zweite Fall vorliegt und somit der ursprüngliche KS-Test nicht verwendet werden kann/darf. Es gibt jedoch entsprechende Anpassungen, welche im Buch von D'Agostino und Stephens beschrieben sind. Da der KS-Test jedoch die geringste Güte von allen EDF-Tests aufweist (vgl. Büning und Trenkler, S. 84), scheint er weniger empfehlenswert.

Btw: Die zweiseitige Variante sollte niemals verwendet werden, da sie verfälscht ist, d. h. die richtige Hypothese kann mit höherer Wahrscheinlichkeit abgelehnt werden als eine falsche ("worser than useless", vgl. S. 196 und S. 23 in Rüger, B. (2002). Statistische Tests, Band 2 aus Test- und Schätztheorie. München, Oldenbourg Wissenschaftsverlag).

Nimm lieber den Anderson-Darling-Test in der Variante für zusammengesetzte Hypothesen oder - noch besser - einen Normalitätstest, d. h. einen Test, welcher speziell für den Fall einer Normalverteilung ausgelegt ist. (Aber nicht den Lilliefors-Test, der aus dem KS-Test hervorgegangen ist.)

Wobei die Grafiken deiner KDE im verlinkten Thread weniger nach Normalverteilung aussahen. Ein Q-Q Plot ist vielleicht die schnellere und einfachere Alternative, wenn keine Signifikanz nachgewiesen werden muss.


Lucky_7
Aktiv
Dabei seit: 21.01.2018
Mitteilungen: 143
Aus:
Beitrag No.5, vom Themenstarter, eingetragen 2019-07-11 22:03

Ich habe festgestellt, dass die Frage doch etwas verwirrend gestellt war.
Ich versuche es noch einmal:

Ich habe insgesamt 300 Datenpunkte gegeben. Durch eine Kerndichteschätzung (KDE) ist es mir gelungen, die Wahrscheinlichkeitsdichtefunktion (PDF) der Verteilung zu schätzen, aus der die 300 Datenpunkte "gezogen" wurden. Das ganze sieht dann so aus:


Es könnte also sein, dass die Verteilung, aus der die 300 Datenpunkte stammen, eine Bimodalverteilung ist. Diese Hypothese möchte ich gerne testen. Bei der Gelegenheit teste ich auch gleich mit, ob die Datenpunkte einer Normalverteilung entspringen.

Im einzelnen teste ich folgendes:
Ist die geschätzte pdf normalverteilt:
python
stat_gauss, pvalue_gauss = kstest(pdf, 'norm', args=(mean(pdf), np.std(pdf)))
, wobei pdf natürlich die geschätzte PDF aus den Datenpunkten entspricht.
Und Ich erhalte als D-Wert: 0.324

Weiterhin teste ich:
python
stat_bimodal, pvalue_bimodal = kstest(pdf, cdf=bimodal_cdf, args=tuple(popt_bimodal))
woebi bimodal_cdf eine selbst geschriebene Funktion beschreibt:
python
def bimodal_cdf(x, weight1, mu1, sigma1, mu2, sigma2): #in a bimodal distribution the weights should add up to 1 --> weight2 = 1-weight1
    """
    CDF of a mixture of two normal distributions.
    """
    return (weight1*norm.cdf(x, mu1, sigma1) +
            (1 - weight1)*norm.cdf(x, mu2, sigma2))
und popt_bimodal sind die Parameter für diese Funktion, die durch curve_fit ermittelt wurden.

Hierbei erhalte ich einen D-Wert von: 0.363

Allerdings macht das alles gar keinen Sinn, weil die entsprechenden CDFs folgendermaßen aussehen:



Der D-Wert für die Bimodalverteilung sollte also kleiner sein. Und der D-Wert für die Gaußverteilung sollte größer sein.

Ist irgendetwas also schief gelaufen?


b_p
Senior
Dabei seit: 11.09.2012
Mitteilungen: 541
Aus: Köln
Beitrag No.6, eingetragen 2019-07-12 20:56

2019-07-11 22:03 - Lucky_7 in Beitrag No. 5 schreibt:
python
stat_gauss, pvalue_gauss = kstest(pdf, 'norm', args=(mean(pdf), np.std(pdf)))
, wobei pdf natürlich die geschätzte PDF aus den Datenpunkten entspricht.

Was ist pdf genau? Sofern es sich nicht um die Liste der Realisationen handeln sollte, hätten wir den Fehler bereits gefunden.

Ungeachtet dessen, machst du aber genau den Fehler, den ich in meinem vorherigen Beitrag beschrieben habe: (zweiseitiger) KS-Test bei zusammengesetzter Hypothese unter Beibehaltung der für diesen Fall nicht vorgesehenen kritischen Werte, d. h. das Ergebnis wird nicht verwertbar sein.




Dieses Forumbeitrag kommt von Matroids Matheplanet
https://https://matheplanet.de

Die URL für dieses Forum-Thema ist:
https://https://matheplanet.de/default3.html?topic=241415=5101
Druckdatum: 2019-08-24 14:06