Wenn wir lineare Regressionsmodelle erstellen, möchten wir manchmal die Konfidenzintervalle für unsere Achsenabschnitte und Koeffizienten kennen. Dies kann folgende und weitere Gründe haben:
- Schätzung: Geben Sie eine Reihe plausibler Werte für die Koeffizienten an, um die Unsicherheit der Schätzungen zu quantifizieren.
- Hypothesentest: Intervalle ungleich Null weisen auf eine signifikante Beziehung zwischen Variablen hin.
- Modellbewertung: Die Intervallbreite spiegelt die Präzision der Koeffizienten wider und hilft bei der Modellauswahl.
- Vergleich: Überlappende Intervalle deuten auf ähnliche Effekte hin, während nicht überlappende Intervalle auf signifikante Unterschiede zwischen den Prädiktoren hinweisen.
Scikit-Study enthält jedoch standardmäßig keine Konfidenzintervalle Statistikmodelle tut. Für meine eigene Regressionsmodellimplementierung mit Numba werde ich daher den Statsmodel-Ansatz als Maßstab verwenden.
Zumindest für mich ist es nicht ganz einfach, der Statistikmodellbibliothek zu folgen, additionally habe ich sie durchgesehen und rückwärts gearbeitet und hoffentlich die Schlüsselelemente gefunden.
Die Lösung finden
In der Dokumentation heißt es:
Das Konfidenzintervall basiert auf der Scholar-t-Verteilung.
In der Statistik werden Konfidenzintervalle für den Mittelwert (oder eine andere Statistik) wie folgt angegeben:
Mittelwert ± SE * z_alpha
Z_alpha ist der kritische Wert basierend auf dem Konfidenzintervall. Bei Fällen mit kleiner Stichprobengröße (was manchmal bei der Regressionsmodellierung der Fall sein kann) verwenden wir jedoch die Scholar-t-Verteilung als Näherungswert. SE ist der Standardfehler des Mittelwerts, additionally die Abweichung des Stichprobenmittelwerts.
Es gibt keine weiteren Particulars, aber ich bin auf den Quellcode gestoßen Hier und die einzige Ergänzung dazu ist so ziemlich:
ci = tremendous(RegressionResults, self).conf_int(alpha=alpha, cols=cols)
Nicht sehr hilfreich! Aber ich habe die Bibliothek dafür gefunden Regressionsergebnisse Klasse weiter GitHub was uns nun ergibt:
decrease = params - q * bse
higher = params + q * bse
Dies ähnelt allmählich der Definition, die wir zuvor gegeben haben – additionally tauchen wir etwas tiefer ein. Wir können das sehen Q stellt den t-Wert des Schülers dar, der sich ergibt aus:
dist = stats.t
q = dist.ppf(1 - alpha / 2, df_resid)
Wobei stats.t das ist SciPy-Paket für die Zufallsvariable Scholar t. Alpha bezieht sich auf das Signifikanzniveau, sodass ein Wert von 0,05 einem zweiseitigen Signifikanzniveau von 95 % entspricht. Ich bin mir nicht sicher, was df_resid Deshalb gehe ich zurück zur Statsmodel-Web site und sehe, dass sie sich auf Folgendes bezieht:
Der verbleibende Freiheitsgrad.
Der Freiheitsgrad ist definiert als die Anzahl der Beobachtungen minus dem Rang der Regressormatrix.
Hier bezieht sich die Anzahl der Beobachtungen auf Datenzeilen im Trainingssatz. Für den Rang der Matrix habe ich gefunden Das hilfreiche Definition:
Warum den Rang ermitteln?
Der Rang verrät uns viel über die Matrix.
Es ist nützlich, um uns mitzuteilen, ob wir eine Probability haben, ein System linearer Gleichungen zu lösen: Wenn der Rang gleich der Anzahl der Variablen ist, können wir möglicherweise eine eindeutige Lösung finden.
Additionally im Wesentlichen Rang sagt uns, wie viele Options „linear unabhängig“ sind und nicht nur Transformationen anderer Options (z. B. ist ein Characteristic doppelt so groß wie ein anderes). In diesem Fall gilt für die meisten wohlgeformten Regressionsdaten: Rang entspricht einfach der Anzahl der Options.
Wir wissen, dass die Scholar-t-Verteilung umso mehr der Normalverteilung ähnelt, je höher die Freiheitsgrade sind. Daher ist es sinnvoll, dass mehr Beobachtungen und weniger Merkmale zu höheren Freiheitsgraden führen.
Jetzt sollten wir herausfinden, was „bse“ ist. In der Dokumentation heißt es:
Die Standardfehler der Parameterschätzungen.
Dies wird in der Bibliothek wie folgt berechnet.
np.sqrt(np.diag(self.cov_params()))
Hier cov_params Ist:
Die Varianz/Kovarianz-Matrix kann einen linearen Kontrast haben
geschätzte Parameter oder alle Parameter multipliziert mit der Skalierung, was der Fall sein wird
handelt es sich normalerweise um eine Schätzung von Sigma². Es wird davon ausgegangen, dass es sich bei der Skala um einen Skalar handelt.
In einem Varianz/Kovarianz-Matrix; Die Diagonalwerte (Cjj) repräsentieren die Varianz der j-ten Elemente und der Relaxation (Cij) repräsentiert die Kovarianz zwischen dem i-ten und dem j-ten Aspect. Hier verwendet np.diag nur die Diagonalelemente (Cjj), daher sind wir für die Berechnung des Standardfehlers nur an der Varianz jedes Merkmals und nicht an deren Kovarianz mit anderen Merkmalen interessiert.
Was ist Maßstab?:
Die Skala ist typischerweise der mittlere quadratische Fehler des geschätzten Modells (Sigma²).
Wenn wir uns den Code jedoch etwas genauer ansehen, können wir erkennen, dass der mittlere quadratische Restfehler verwendet wird. Dieser wird wie der mittlere quadratische Fehler berechnet, jedoch:
Die Summe der quadrierten Residuen dividiert durch die Restfreiheitsgrade.
Lösung
Jetzt haben wir additionally alle Teile des Puzzles. Fassen wir es zusammen:
from scipy import stats
import numpy as npalpha # the importance stage
X # regression matrix with intercept, measurement (n_samples, n_features)
y # goal vector, measurement n_samples
y_fitted # fitted values
params # fitted intercept + coefficient weights
n_samples = X.form[0]
rank = np.linalg.matrix_rank(X)
df_resid = n_samples - rank # levels of freedom
residuals = y_fitted - y
scale = np.sum(residuals ** 2) / df_resid # resid mse
C = np.linalg.inv(X.T @ X) * scale # covar matrix
v = np.diag(C) # parameter variances
bse = np.sqrt(v) # customary errors
q = stats.t.ppf(1 - alpha / 2, df_resid) # t-value
decrease = params - q * bse
higher = params + q * bse
Ergebnisse
Lassen Sie uns nun diese Implementierung im Vergleich zu Statistikmodellen überprüfen:
from sklearn.linear_model import LinearRegression
from sklearn.datasets import make_regression
import statsmodels.api as sm
from scipy import stats
import pandas as pd
import numpy as np# create random knowledge
X, y = make_regression(n_samples=150, n_features=15, noise=10)
significance_level = .05
# match statsmodels
sm_model = sm.OLS(y, sm.add_constant(X))
outcomes = sm_model.match()
sm_ci = outcomes.conf_int(alpha=significance_level)
# match scikit-learn
sk_model = LinearRegression(fit_intercept=True)
sk_model.match(X, y)
y_fitted = sk_model.predict(X)
# add intercept
intercept = np.ones(form=(X.form[0], 1), dtype=np.float64)
X = np.concatenate([intercept, X], axis=1)
# residuals levels of freedom
n_samples = X.form[0]
rank = np.linalg.matrix_rank(X)
df_resid = n_samples - rank
# t-value
q = stats.t.ppf(1 - significance_level / 2, df_resid)
# scale
residuals = y_fitted - y
scale = np.sum(residuals ** 2) / df_resid
# customary errors
C = np.linalg.inv(X.T @ X) * scale
v = np.diag(C)
bse = np.sqrt(v)
# prepend intercept to coeffs
params = np.r_[m.intercept_, m.coef_]
labels = ['intercept'] + checklist(map(lambda c: 'coef_' + str(c), vary(1, X.form[1])))
outcomes = pd.DataFrame(
{'label': labels,
'lower_bound': params - q * bse,
'params': params + q * bse,
'upper_bound': higher,
'lower_statsmodels': sm_ci[:, 0],
'upper_statsmodels': sm_ci[:, 1]
}
)
Hyperlink zu GitHub: jcatankard/NumbaML: Regressionsmodellierung in Numba implementiert
