Analyse einer Stichprobe#
Zuerst erstellen wir einige Zufallsvariablen. Wir setzen einen Seed, damit wir bei jedem Durchlauf identische Ergebnisse zum Betrachten erhalten. Als Beispiel nehmen wir eine Stichprobe aus der Student'schen t-Verteilung.
>>> import numpy as np
>>> import scipy.stats as stats
>>> x = stats.t.rvs(10, size=1000)
Hier setzen wir den erforderlichen Formparameter der t-Verteilung, der in der Statistik den Freiheitsgraden entspricht, auf 10. Mit size=1000 bedeutet, dass unsere Stichprobe aus 1000 unabhängig gezogenen (Pseudo-)Zufallszahlen besteht. Da wir die Schlüsselwortargumente loc und scale nicht spezifiziert haben, sind diese auf ihre Standardwerte Null und Eins gesetzt.
Deskriptive Statistiken#
x ist ein NumPy-Array und wir haben direkten Zugriff auf alle Array-Methoden, z.B.:
>>> print(x.min()) # equivalent to np.min(x)
-3.78975572422 # random
>>> print(x.max()) # equivalent to np.max(x)
5.26327732981 # random
>>> print(x.mean()) # equivalent to np.mean(x)
0.0140610663985 # random
>>> print(x.var()) # equivalent to np.var(x))
1.28899386208 # random
Wie vergleichen sich die Stichprobeneigenschaften mit ihren theoretischen Gegenstücken?
>>> m, v, s, k = stats.t.stats(10, moments='mvsk')
>>> n, (smin, smax), sm, sv, ss, sk = stats.describe(x)
>>> sstr = '%-14s mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f'
>>> print(sstr % ('distribution:', m, v, s ,k))
distribution: mean = 0.0000, variance = 1.2500, skew = 0.0000, kurtosis = 1.0000 # random
>>> print(sstr % ('sample:', sm, sv, ss, sk))
sample: mean = 0.0141, variance = 1.2903, skew = 0.2165, kurtosis = 1.0556 # random
Hinweis: stats.describe verwendet den unverzerrten Schätzer für die Varianz, während np.var der verzerrte Schätzer ist.
Für unsere Stichprobe weichen die Stichprobenstatistiken geringfügig von ihren theoretischen Gegenstücken ab.
T-Test und KS-Test#
Wir können den t-Test verwenden, um zu prüfen, ob sich der Mittelwert unserer Stichprobe statistisch signifikant von der theoretischen Erwartung unterscheidet.
>>> print('t-statistic = %6.3f pvalue = %6.4f' % stats.ttest_1samp(x, m))
t-statistic = 0.391 pvalue = 0.6955 # random
Der p-Wert beträgt 0,7, das bedeutet, dass wir bei einem Alpha-Fehler von z.B. 10 % die Hypothese, dass der Stichprobenmittelwert gleich Null ist, der Erwartung der Standard-t-Verteilung, nicht ablehnen können.
Als Übung können wir unseren t-Test auch direkt berechnen, ohne die bereitgestellte Funktion zu verwenden, was uns dasselbe Ergebnis liefern sollte, und das tut es auch.
>>> tt = (sm-m)/np.sqrt(sv/float(n)) # t-statistic for mean
>>> pval = stats.t.sf(np.abs(tt), n-1)*2 # two-sided pvalue = Prob(abs(t)>tt)
>>> print('t-statistic = %6.3f pvalue = %6.4f' % (tt, pval))
t-statistic = 0.391 pvalue = 0.6955 # random
Der Kolmogorov-Smirnov-Test kann verwendet werden, um die Hypothese zu testen, dass die Stichprobe aus der Standard-t-Verteilung stammt.
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 't', (10,)))
KS-statistic D = 0.016 pvalue = 0.9571 # random
Auch hier ist der p-Wert hoch genug, dass wir die Hypothese, dass die Zufallsstichprobe tatsächlich nach der t-Verteilung verteilt ist, nicht ablehnen können. In realen Anwendungen wissen wir nicht, wie die zugrundeliegende Verteilung ist. Wenn wir den Kolmogorov-Smirnov-Test unserer Stichprobe gegen die Standard-Normalverteilung durchführen, können wir ebenfalls die Hypothese nicht ablehnen, dass unsere Stichprobe durch die Normalverteilung generiert wurde, da in diesem Beispiel der p-Wert fast 40 % beträgt.
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 'norm'))
KS-statistic D = 0.028 pvalue = 0.3918 # random
Die Standard-Normalverteilung hat jedoch eine Varianz von 1, während unsere Stichprobe eine Varianz von 1,29 hat. Wenn wir unsere Stichprobe standardisieren und gegen die Normalverteilung testen, ist der p-Wert wieder groß genug, dass wir die Hypothese, dass die Stichprobe aus der Normalverteilung stammt, nicht ablehnen können.
>>> d, pval = stats.kstest((x-x.mean())/x.std(), 'norm')
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % (d, pval))
KS-statistic D = 0.032 pvalue = 0.2397 # random
Hinweis: Der Kolmogorov-Smirnov-Test geht davon aus, dass wir gegen eine Verteilung mit gegebenen Parametern testen. Da im letzten Fall Mittelwert und Varianz geschätzt wurden, ist diese Annahme verletzt und die Verteilung der Teststatistik, auf der der p-Wert basiert, ist nicht korrekt.
Schwänze der Verteilung#
Schließlich können wir den oberen Schwanz der Verteilung überprüfen. Wir können die Perzentilpunktfunktion ppf, die die Umkehrung der CDF-Funktion ist, verwenden, um die kritischen Werte zu erhalten, oder direkter die Umkehrung der Überlebensfunktion.
>>> crit01, crit05, crit10 = stats.t.ppf([1-0.01, 1-0.05, 1-0.10], 10)
>>> print('critical values from ppf at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' % (crit01, crit05, crit10))
critical values from ppf at 1%, 5% and 10% 2.7638 1.8125 1.3722
>>> print('critical values from isf at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' % tuple(stats.t.isf([0.01,0.05,0.10],10)))
critical values from isf at 1%, 5% and 10% 2.7638 1.8125 1.3722
>>> freq01 = np.sum(x>crit01) / float(n) * 100
>>> freq05 = np.sum(x>crit05) / float(n) * 100
>>> freq10 = np.sum(x>crit10) / float(n) * 100
>>> print('sample %%-frequency at 1%%, 5%% and 10%% tail %8.4f %8.4f %8.4f' % (freq01, freq05, freq10))
sample %-frequency at 1%, 5% and 10% tail 1.4000 5.8000 10.5000 # random
In allen drei Fällen hat unsere Stichprobe mehr Gewicht im oberen Schwanz als die zugrundeliegende Verteilung. Wir können kurz eine größere Stichprobe überprüfen, um zu sehen, ob wir eine engere Übereinstimmung erhalten. In diesem Fall ist die empirische Häufigkeit ziemlich nah an der theoretischen Wahrscheinlichkeit, aber wenn wir dies mehrmals wiederholen, sind die Schwankungen immer noch ziemlich groß.
>>> freq05l = np.sum(stats.t.rvs(10, size=10000) > crit05) / 10000.0 * 100
>>> print('larger sample %%-frequency at 5%% tail %8.4f' % freq05l)
larger sample %-frequency at 5% tail 4.8000 # random
Wir können sie auch mit dem Schwanz der Normalverteilung vergleichen, die weniger Gewicht in den Schwänzen hat.
>>> print('tail prob. of normal at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' %
... tuple(stats.norm.sf([crit01, crit05, crit10])*100))
tail prob. of normal at 1%, 5% and 10% 0.2857 3.4957 8.5003
Der Chi-Quadrat-Test kann verwendet werden, um zu prüfen, ob bei einer endlichen Anzahl von Bins die beobachteten Häufigkeiten signifikant von den Wahrscheinlichkeiten der hypothesierten Verteilung abweichen.
>>> quantiles = [0.0, 0.01, 0.05, 0.1, 1-0.10, 1-0.05, 1-0.01, 1.0]
>>> crit = stats.t.ppf(quantiles, 10)
>>> crit
array([ -inf, -2.76376946, -1.81246112, -1.37218364, 1.37218364,
1.81246112, 2.76376946, inf])
>>> n_sample = x.size
>>> freqcount = np.histogram(x, bins=crit)[0]
>>> tprob = np.diff(quantiles)
>>> nprob = np.diff(stats.norm.cdf(crit))
>>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
>>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
>>> print('chisquare for t: chi2 = %6.2f pvalue = %6.4f' % (tch, tpval))
chisquare for t: chi2 = 2.30 pvalue = 0.8901 # random
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 64.60 pvalue = 0.0000 # random
Wir sehen, dass die Standard-Normalverteilung klar abgelehnt wird, während die Standard-t-Verteilung nicht abgelehnt werden kann. Da die Varianz unserer Stichprobe von beiden Standardverteilungen abweicht, können wir den Test erneut durchführen und die Schätzung für Skala und Lage berücksichtigen.
Die Fit-Methode der Verteilungen kann verwendet werden, um die Parameter der Verteilung zu schätzen, und der Test wird mit Wahrscheinlichkeiten der geschätzten Verteilung wiederholt.
>>> tdof, tloc, tscale = stats.t.fit(x)
>>> nloc, nscale = stats.norm.fit(x)
>>> tprob = np.diff(stats.t.cdf(crit, tdof, loc=tloc, scale=tscale))
>>> nprob = np.diff(stats.norm.cdf(crit, loc=nloc, scale=nscale))
>>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
>>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
>>> print('chisquare for t: chi2 = %6.2f pvalue = %6.4f' % (tch, tpval))
chisquare for t: chi2 = 1.58 pvalue = 0.9542 # random
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 11.08 pvalue = 0.0858 # random
Unter Berücksichtigung der geschätzten Parameter können wir die Hypothese, dass unsere Stichprobe aus einer Normalverteilung stammt (auf dem 5 %-Niveau), immer noch ablehnen, aber wieder, mit einem p-Wert von 0,95, können wir die t-Verteilung nicht ablehnen.
Spezielle Tests für Normalverteilungen#
Da die Normalverteilung die häufigste Verteilung in der Statistik ist, gibt es mehrere zusätzliche Funktionen, um zu testen, ob eine Stichprobe aus einer Normalverteilung gezogen worden sein könnte.
Zuerst können wir testen, ob Schiefe und Kurtosis unserer Stichprobe signifikant von denen einer Normalverteilung abweichen.
>>> print('normal skewtest teststat = %6.3f pvalue = %6.4f' % stats.skewtest(x))
normal skewtest teststat = 2.785 pvalue = 0.0054 # random
>>> print('normal kurtosistest teststat = %6.3f pvalue = %6.4f' % stats.kurtosistest(x))
normal kurtosistest teststat = 4.757 pvalue = 0.0000 # random
Diese beiden Tests werden im Normalitätstest kombiniert.
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' % stats.normaltest(x))
normaltest teststat = 30.379 pvalue = 0.0000 # random
In allen drei Tests sind die p-Werte sehr niedrig und wir können die Hypothese ablehnen, dass unsere Stichprobe die Schiefe und Kurtosis der Normalverteilung aufweist.
Da Schiefe und Kurtosis unserer Stichprobe auf zentralen Momenten basieren, erhalten wir genau dieselben Ergebnisse, wenn wir die standardisierte Stichprobe testen.
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest((x-x.mean())/x.std()))
normaltest teststat = 30.379 pvalue = 0.0000 # random
Da die Normalität so stark abgelehnt wird, können wir prüfen, ob der Normaltest für andere Fälle vernünftige Ergebnisse liefert.
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest(stats.t.rvs(10, size=100)))
normaltest teststat = 4.698 pvalue = 0.0955 # random
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest(stats.norm.rvs(size=1000)))
normaltest teststat = 0.613 pvalue = 0.7361 # random
Beim Test auf Normalität einer kleinen Stichprobe von t-verteilten Beobachtungen und einer großen Stichprobe von normalverteilten Beobachtungen können wir in keinem der Fälle die Nullhypothese ablehnen, dass die Stichprobe aus einer Normalverteilung stammt. Im ersten Fall liegt dies daran, dass der Test nicht mächtig genug ist, um eine t- und eine normalverteilte Zufallsvariable in einer kleinen Stichprobe zu unterscheiden.