Wenn Sie eine Datei wie Python programmgesteuert lesen oder verarbeiten, sollten Sie sich angewöhnen, zu prüfen, ob die Datei wie erwartet gelesen werden kann. Wenn Sie im Fall von Jupyter Notebook nur den Variablennamen eingeben, wird der Inhalt dieser Variablen übrigens in Out angezeigt, was praktisch ist.

chem

79 rows × 9 columns

sample

79 rows × 66 columns

Kapitel 3 Datenaufbereitung

Nach dem Lesen der Daten besteht der nächste Schritt darin, die Daten vorzubereiten.
Zunächst werden die zweigeteilten Daten (chem und Probe) integriert und nur die erforderlichen Daten extrahiert. Dieses Mal möchten wir die Daten von Katsuos ΣPCB verwenden, daher werden nur die Daten extrahiert, deren Wert in der Spalte "ChemicalName" "ΣPCB" ist.

df = pd.merge(chem, sample, on="SampleID")
data = df[df["ChemicalName"] == "ΣPCBs"]

Überprüfen Sie als Nächstes, ob die Maßeinheit unterschiedlich ist. Dies liegt daran, dass es nicht möglich ist, einfach zu vergleichen und zu integrieren, wenn die Dateneinheiten wie in Teil 1 unterschiedlich sind.

data["Unit"].unique()

array(['ng/g wet'], dtype=object)

​

Sie können die eindeutige Methode

pandas verwenden, um eine Liste der in diesem Datenrahmen enthaltenen Werte anzuzeigen. Wenn Sie hier eine Liste von Werten ausgeben, die in der Spalte "Einheit" enthalten sind, sehen Sie, dass es sich nur um "ng / g nass" handelt. Diesmal ist es nicht erforderlich, die Daten nach Einheiten zu teilen. </ p>

Schließlich wird die Spalte mit nur N / A gelöscht und die Datenvorbereitung abgeschlossen.

data = data.dropna(how='all', axis=1)
data

79 rows × 35 columns

Kap.4 Punktschätzung

Von den statistischen Schätzungen, die die Population aus der Stichprobe schätzen, ist die Punktschätzung die punktgenaue Schätzung des Werts. Hier wird die ΣPCB-Konzentration der Population (ganzes Individuum im Sammelgebiet) aus der Stichprobe der im Bonito nachgewiesenen ΣPCB-Konzentration geschätzt.

Überprüfen Sie zunächst, wie viele Jahre Daten enthalten sind, um jedes Jahr zu berechnen, und sehen Sie den Übergang der Änderung.

data['CollectionYear'].unique()

array([1998, 1997, 1999, 2001], dtype=int64)

Aus der obigen eindeutigen Methode wurde festgestellt, dass der Datensatz Daten für die drei Jahre 1997-1999 und 2001 enthält. Lassen Sie uns zunächst die Daten von 1997 herausnehmen. Ändern Sie zu diesem Zeitpunkt die extrahierten Daten in das Numpy-Format ndarray ^{1 , damit zukünftige Berechnungen einfacher werden.}

pcb_1997 = np.array(data[data['CollectionYear']==1997]["MeasuredValue"]) #Extrahieren Sie nur 1997 Messungen
pcb_1997

array([ 10.72603788,   9.22208078,   7.59790835,  30.95079465,
        15.27462553,  14.15719633,  13.28955903,  14.87712806,
         9.86650189,  18.26554514,   3.39951845,   6.58172781,
        12.43564814,   6.1948639 ,   6.41605666,   4.98827291,
        12.36669815,  31.17955551,   8.16184346,   4.60893266,
        36.85826409,  52.99841724,  39.22500351,  53.92302702,
        69.4308048 ,  73.97686479, 125.3887794 ,  45.39974771,
        54.12726127,  39.77794045, 101.2736126 ,  38.06220403,
       126.8301693 ,  70.25308435,  31.24246301,  21.3958656 ,
        41.85726522,  30.91112132,  81.12597135,  10.76755148,
        24.20442213,  24.57497594,  14.84353549,  59.53687389,
        52.78443082,   8.4644697 ,   4.15293758,   3.31957452,
         4.51832675,   6.98373973])

In ähnlicher Weise werden auch die Daten für 1998 und 1999 extrahiert.

pcb_1998 = np.array(data[data['CollectionYear']==1998]["MeasuredValue"]) #Extrahieren Sie nur Messungen von 1998
pcb_1999 = np.array(data[data['CollectionYear']==1999]["MeasuredValue"]) #Extrahieren Sie nur Messungen von 1999

Hier werden die unverzerrten Schätzungen von Mittelwert und Varianz berechnet.
Erstens die unvoreingenommene Schätzung des Mittelwerts ($ \ hat {\ mu} $), die die Tatsache ausnutzt, dass der erwartete Wert des Stichprobenmittelwerts ($ \ override {X} $) gleich dem Populationsmittelwert ist. (Siehe Formel unten)
$$ E \left(\overline X \right) = E\left(\frac{1}{n} \sum_{i=1}^{n} \left(x_i\right)\right) = \frac{1}{n}\sum_{i=1}^{n} E\left(x_i\right) = \frac{1}{n} \times n\mu = \mu \\ \therefore \hat\mu = \overline{x} $$

s_mean_1997 = np.mean(pcb_1997)
s_mean_1997

31.775384007760003

Berechnen Sie auf ähnliche Weise die unverzerrte Schätzung der Varianz.
Zu diesem Zeitpunkt nimmt der erwartete Wert der Stichprobenvarianz ($ S ^ 2 $) nicht den gleichen Wert wie die Populationsvarianz ($ \ sigma ^ 2 $) an, und es ist erforderlich, stattdessen die unverzerrte Varianz ($ s ^ 2 $) zu erhalten. Beachten.
$$\hat\sigma^2 \neq S^2 = \frac{1}{n} \sum_{i=1}^{n} \left( x_i - \overline X \right) \\ \hat\sigma^2 = s^2 = \frac{1}{n-1} \sum_{i=1}^{n} \left( x_i - \overline X \right)$$

Beachten Sie, dass die var-Funktion von Numpy beide Varianzen berechnen kann und die unverzerrte Varianz mit dem Parameter ddof = 1 ausgegeben wird. Beachten Sie jedoch, dass die Stichprobenverteilung von ddof = 0 standardmäßig ausgegeben wird.
$$\mathrm{np.var}\left(x_1 \ldots x_n, \mathrm{ddof=0}\right): S^2 = \frac{1}{n} \sum_{i=1}^{n} \left( x_i - \overline X \right) \\ \mathrm{np.var}\left(x_1 \ldots x_n, \mathrm{ddof=1}\right): \hat\sigma^2 = \frac{1}{n-1} \sum_{i=1}^{n} \left( x_i - \overline X \right)$$

u_var_1997 = np.var(pcb_1997, ddof=1)
u_var_1997 

942.8421749786518

Berechnen Sie in ähnlicher Weise die unvoreingenommenen Schätzungen des Mittelwerts und der Varianz für 1998 und 1999.

s_mean_1998, s_mean_1999 = np.mean(pcb_1998), np.mean(pcb_1999)
u_var_1998, u_var_1999 = np.var(pcb_1998, ddof=1), np.var(pcb_1999, ddof=1)

Hier fassen wir die erhaltenen repräsentativen Werte zusammen. Erstens

s_mean_1997, s_mean_1998, s_mean_1999

(31.775384007760003, 17.493267312533337, 30.583242522000003)

u_var_1997, u_var_1998, u_var_1999

(942.8421749786518, 240.2211176248311, 1386.7753819003349)

Kapitel 5 Abschnitt Schätzung und Vertrauensabschnitt

In Kapitel 5 wird im Gegensatz zu der in Kapitel 4 erhaltenen Punktschätzung die Intervallschätzung durchgeführt, um den Populationsmittelwert und die Populationsvarianz innerhalb eines statistisch konstanten Bereichs zu schätzen.

Abschnitt 5-1 Intervallschätzung des Bevölkerungsdurchschnitts

Überprüfen Sie zunächst die Anzahl der Daten im Datensatz jedes Jahres, bevor Sie das Intervall schätzen. Die Anzahl der Daten kann mithilfe der in Python standardmäßig implementierten len-Funktion berechnet werden.

n_1997 = len(pcb_1997)
n_1997

50

n_1998, n_1999 = len(pcb_1998), len(pcb_1999)
n_1998, n_1999

(15, 13)

Aus dem oben Gesagten ist die Gesamtzahl der Daten im Datensatz für jedes Jahr von 1997 bis 1999 bekannt.

Von diesen ist der Datensatz von 1997 eine große Stichprobe von $ n = 50 $, und die Datensätze von 1998 und 1999 sind $ n = \ left \ {\ begin {array} {ll}. 15 & \ left (1998 \ right) \\ 13 & \ left (1999 \ right) \ end {array} \ right. $, Eine kleine Stichprobe.
Daher ist zu beachten, dass die Verarbeitung der Intervallschätzung danach geringfügig anders ist.

Zunächst wird das Intervall des Bevölkerungsdurchschnitts aus dem Datensatz von 1997 geschätzt. In diesem Fall ist der Stichprobenmittelwert ($ \ override X $) eine Normalverteilung $ N \ left (\ mu, \ frac {s ^ 2} {, da die Populationsvarianz unbekannt ist und eine große Stichprobe ($ n> 30 $) n} \ right) $ kann angenähert werden. Wenn daher der Populationsmittelwert durch die Zuverlässigkeit ($ \ alpha $) geschätzt wird, ist das Konfidenzintervall wie in der folgenden Gleichung gezeigt. $\overline X - z_\frac{\alpha}{2} \sqrt{\frac{s^2}{n}} < \mu < \overline X - z_\frac{\alpha}{2} \sqrt{\frac{s^2}{n}} $

Wenn Sie in Python stas.norm.interval () von Scipy verwenden, können Sie den Bereich von Alpha × 100% in der Normalverteilung von Mittelwert (loc) und Standardabweichung (Skala) erhalten, die auf dem Median zentriert ist.
Hier wird das Konfidenzintervall durch die Zuverlässigkeit berechnet ($ \ alpha = 0,95 $).

m_interval_1997 = stats.norm.interval(alpha=0.95, loc=s_mean_1997, scale=math.sqrt(pcb_1997.var(ddof=1)/n_1997))
m_interval_1997

(23.26434483549182, 40.28642318002819)

Als Nächstes schätzen wir das mittlere Populationsintervall für die Datensätze 1998 und 1999. Dies sind kleine Stichproben ($ n \ leq 30 $) mit unbekannter Populationsvarianz. In diesem Fall ist der Populationsmittelwert $ \ mu $ nicht die Normalverteilung $ N \ left (\ mu, \ frac {s ^ 2} {n} \ right) $, sondern die t-Verteilung mit Freiheitsgraden ($ n-1 $). Benutzen. Wenn daher der Populationsmittelwert durch die Zuverlässigkeit ($ \ alpha $) geschätzt wird, ist das Konfidenzintervall wie in der folgenden Gleichung gezeigt. $\overline X - t_\frac{\alpha}{2}\left(n-1\right)\sqrt{\frac{s^2}{n}} < \mu < \overline X + t_\frac{\alpha}{2}\left(n-1\right)\sqrt{\frac{s^2}{n}} $

Wenn Sie in Python stats.t.interval () verwenden, wird der Bereich, in dem die t-Verteilung von Mittelwert (loc), Standardabweichung (Skala) und Freiheitsgrad (df) Alpha × 100% beträgt, auf den Medianwert zentriert. Erhältlich als.
Hier wird das Konfidenzintervall basierend auf der Zuverlässigkeit berechnet ($ \ alpha = 0,95 $).

m_interval_1998 = stats.t.interval(alpha=0.95, df=n_1998-1, loc=s_mean_1998, scale=math.sqrt(pcb_1998.var(ddof=1)/n_1998))
m_interval_1999 = stats.t.interval(alpha=0.95, df=n_1999-1, loc=s_mean_1999, scale=math.sqrt(pcb_1999.var(ddof=1)/n_1999))

m_interval_1997, m_interval_1998, m_interval_1999

((23.26434483549182, 40.28642318002819),
 (8.910169386248537, 26.076365238818138),
 (8.079678286109523, 53.086806757890486))

Das 95% -Konfidenzintervall bedeutet, dass der Bevölkerungsdurchschnitt mit einer Wahrscheinlichkeit von 95% innerhalb dieses Bereichs liegt. Das heißt, wenn die Zuverlässigkeit ($ \ alpha $) verringert wird, wird die Wahrscheinlichkeit verringert, dass der Bevölkerungsdurchschnitt in das Konfidenzintervall aufgenommen wird, und gleichzeitig wird das Konfidenzintervall verringert.

stats.norm.interval(alpha=0.9, loc=s_mean_1997, scale=math.sqrt(pcb_1997.var(ddof=1)/n_1997))

(24.63269477364296, 38.91807324187704)

Abschnitt 5-2 Abschnittsschätzung der Populationsvarianz

Als nächstes wird die Intervallschätzung der Populationsvarianz durchgeführt. Bei der Intervallschätzung der Populationsvarianz ($ \ sigma ^ 2 $) hat $ \ frac {\ left (n-1 \ right) s ^ 2} {\ sigma ^ 2} $ einen Freiheitsgrad $ (n-1) $. Verwenden Sie diese Option, um der Verteilung von $ \ chi ^ 2 $ von zu folgen. $\chi_\frac{\alpha}{2}\left(n-1\right) \leq \frac{\left( n-1 \right)s^2}{\sigma^2} \leq \chi_{1-\frac{\alpha}{2}}\left(n-1\right)$

Zuerst der Prozentpunkt ($ \ chi_ \ frac {\ alpha} {2} \ left (n-1 \ right), \ chi_ der $ \ chi ^ 2 $ -Verteilung der Freiheit ($ n-1 $) Suchen Sie {1- \ frac {\ alpha} {2}} \ left (n-1 \ right) $). Hier beträgt die Zuverlässigkeit 0,95.
In Python können Sie mit Scipys stats.chi2.interval () den Bereich von Alpha x 100% des Freiheitsgrades (df) ermitteln.

chi_025_1997, chi_975_1997 = stats.chi2.interval(alpha=0.95, df=n_1997-1)
chi_025_1997, chi_975_1997

(31.554916462667137, 70.22241356643451)

Suchen Sie als Nächstes das Konfidenzintervall. Beziehen Sie sich für die Ableitung auf die folgende Formel. $\frac{\left(n-1\right)s^2}{\chi_\frac{\alpha}{2}\left(n-1\right)} \leq \sigma^2 \leq \frac{\left(n-1\right)s^2}{\chi_{1-\frac{\alpha}{2}}\left(n-1\right)}$

v_interval_1997 = (n_1997 - 1)*np.var(pcb_1997, ddof=1) / chi_975_1997, (n_1997 - 1)*np.var(pcb_1997, ddof=1) / chi_025_1997
v_interval_1997

(657.8991553778869, 1464.0909168183869)

In ähnlicher Weise wird für die Datensätze 1998 und 1999 die Intervallintervallschätzung der Varianz durchgeführt. Beachten Sie, dass im Gegensatz zur Intervallschätzung des Populationsmittelwerts die Verteilung von $ \ frac {\ left (n-1 \ right) s ^ 2} {\ sigma ^ 2} $ unabhängig von der Stichprobengröße $ \ chi ^ 2 beträgt. Folgen Sie der $ -Verteilung.

chi_025_1998, chi_975_1998 = stats.chi2.interval(alpha=0.95, df=n_1998-1)
chi_025_1999, chi_975_1999 = stats.chi2.interval(alpha=0.95, df=n_1999-1)

v_interval_1998 = (n_1998 - 1)*np.var(pcb_1998, ddof=1) / chi_975_1998, (n_1998 - 1)*np.var(pcb_1998, ddof=1) / chi_025_1998
v_interval_1999 = (n_1999 - 1)*np.var(pcb_1999, ddof=1) / chi_975_1999, (n_1999 - 1)*np.var(pcb_1999, ddof=1) / chi_025_1999

v_interval_1997, v_interval_1998, v_interval_1999

((657.8991553778869, 1464.0909168183869),
 (128.76076176378118, 597.4878836139195),
 (713.0969734866349, 3778.8609867211235))

chi_025_1997, chi_975_1997 = stats.chi2.interval(alpha=0.9, df=n_1997-1)
(n_1997 - 1)*np.var(pcb_1997, ddof=1) / chi_975_1997, (n_1997 - 1)*np.var(pcb_1997, ddof=1) / chi_025_1997

(696.4155490924242, 1361.5929987004467)

Kapitel 6 Visualisierung der Schätzergebnisse

Lassen Sie uns nun den in Kapitel 4 und Kapitel 5 geschätzten Bevölkerungsdurchschnitt in einem Diagramm visualisieren.

Zunächst werden die in Kapitel 4 geschätzten Werte des Bevölkerungsmittelwerts in Zeitreihen zusammengefasst.

x_list = [1997, 1998, 1999]
y_list = [s_mean_1997, s_mean_1998, s_mean_1999]

Als nächstes wird das in Kapitel 5 geschätzte Konfidenzintervall des Bevölkerungsdurchschnitts mit einer Zuverlässigkeit von 95% in Zeitreihen zusammengefasst.

interval_list = []
interval_list.append(m_interval_1997)
interval_list.append(m_interval_1998)
interval_list.append(m_interval_1999)
interval_list

[(23.26434483549182, 40.28642318002819),
 (8.910169386248537, 26.076365238818138),
 (8.079678286109523, 53.086806757890486)]

Das 95% -Konfidenzintervall des Bevölkerungsdurchschnitts kann nicht für die Visualisierung verwendet werden, daher wird die Breite des Konfidenzintervalls berechnet.

interval_list = np.array(interval_list).T[1] - y_list
x_list, y_list, interval_list

([1997, 1998, 1999],
 [31.775384007760003, 17.493267312533337, 30.583242522000003],
 array([ 8.51103917,  8.58309793, 22.50356424]))

Visualisieren Sie abschließend mit matplotlib. Das Vertrauensintervall wird mit einer Fehlerleiste angezeigt. Verwenden Sie die Fehlerbalkenmethode, um die Fehlerleiste mit matplotlib anzuzeigen.
Bei dieser Methode werden der X-Achsenwert (hier Jahr), der Y-Achsenwert (hier der Populationsdurchschnitt der Punktschätzung) und die Länge des Fehlerbalkens (hier die Breite des Konfidenzintervalls) angegeben.

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.errorbar(x=x_list, y=y_list, yerr=interval_list, fmt='o-', capsize=4, ecolor='red')
plt.xticks(x_list)
ax.set_title("Katsuwonus pelamis")
ax.set_ylabel("ΣPCBs [ng/g wet]")
plt.show()

Fußnote

^{1 Ein Datenformat, in dem die Matrix n-ter Ordnung in Numpy gespeichert ist.}