Jak obliczyć przedział prognozowania w GLM (Gamma) / TweedieRegression w Pythonie?
Sprawdziłem wiele źródeł z sieci na temat przeprowadzania interwału prediciton, szczególnie w funkcji GLM. Jedno z podejść dotyczy interwałów przewidywania dla uczenia maszynowegohttps://machinelearningmastery.com/prediction-intervals-for-machine-learning/od Jasona Brownlee. Jednak jego metoda jest ukierunkowana na regresję liniową i w pewnym stopniu może nie być odpowiednia dla GLM (Gamma). Innym podejściem, które znalazłem, jest użycie metody ładowania początkowego do przeprowadzenia przedziału predykcji. Jednak obliczenia były bardzo czasochłonne, a pamięć mojego komputera została zabita podczas uruchamiania funkcji z artykułuhttps://saattrupdan.github.io/2020-03-01-bootstrap-prediction/. Nie mam pojęcia, jak prawidłowo przeprowadzić przedział predykcji w GLM (najprawdopodobniej Gamma) w Pythonie zamiast w R. Znalazłem powiązany pakiet w R, ale nie chcę używać R do prowadzenia interwału. Inną powiązaną informacją, którą znalazłem w sieci, jest Gamma GLM - Derive prediction interval for new x_i: Gamma GLM - Derive prediction interval for new x_i .
Odpowiedzi
Jest to trochę skomplikowane, ale powinno być wykonalne.
Jak mówi ten post, aby uzyskać przedział prognozy, musisz zintegrować niepewność ze współczynnikami. Trudno to zrobić analitycznie, ale zamiast tego możemy to zasymulować. Oto kilka danych dotyczących regresji gamma
N = 100
x = np.random.normal(size = N)
true_beta = np.array([0.3])
eta = 0.8 + x*true_beta
mu = np.exp(eta)
shape = 10
#parameterize gamma in terms of shaope and scale
y = gamma(a=shape, scale=mu/shape).rvs()
Teraz dopasuję regresję gamma do tych danych
X = sm.tools.add_constant(x)
gamma_model = sm.GLM(y, X, family=sm.families.Gamma(link = sm.families.links.log()))
gamma_results = gamma_model.fit()
gamma_results.summary()
Generalized Linear Model Regression Results
Dep. Variable: ,y , No. Observations: , 100
Model: ,GLM , Df Residuals: , 98
Model Family: ,Gamma , Df Model: , 1
Link Function: ,log , Scale: ,0.075594
Method: ,IRLS , Log-Likelihood: , -96.426
Date: ,Mon, 30 Nov 2020, Deviance: , 7.7252
Time: ,22:45:07 , Pearson chi2: , 7.41
No. Iterations: ,7 , ,
Covariance Type:,nonrobust , ,
, coef , std err , z ,P>|z| , [0.025 , 0.975]
const, 0.8172, 0.028, 29.264, 0.000, 0.762, 0.872
x1 , 0.2392, 0.029, 8.333, 0.000, 0.183, 0.296
Tak długo, jak mam wystarczającą ilość danych, możemy dokonać normalnego przybliżenia rozkładu próbkowania współczynników.
Średnią i kowariancję można uzyskać z podsumowania modelu.
beta_samp_mean = gamma_results.params
beta_samp_cov = gamma_results.cov_params()
dispersion = gamma_results.scale
Teraz jest tylko kwestia próbkowania fałszywych danych przy użyciu tych szacunków i kwantyli.
X_pred = np.linspace(-2, 2)
X_pred = sm.tools.add_constant(X_pred)
num_samps = 100_000
possible_coefficients = np.random.multivariate_normal(mean = beta_samp_mean, cov = beta_samp_cov, size = num_samps)
linear_predictions = [X_pred@b for b in possible_coefficients]
y_hyp = gamma(a=1/dispersion, scale = np.exp(linear_predictions)*dispersion).rvs()
# Here is the prediction interval
l, u = np.quantile(y_hyp, q=[0.025, 0.975], axis = 0)
Łatwo jest wtedy wykreślić przedział przewidywania
yhat = gamma_results.predict(X_pred)
fig, ax = plt.subplots(dpi = 120)
plt.plot(X_pred[:,1], yhat, color = 'red', label = 'Estimated')
plt.plot(X_pred[:, 1], np.exp(0.8 + X_pred[:, 1]*true_beta), label = 'Truth')
plt.fill_between(X_pred[:, 1], l, u, color = 'red', alpha = 0.1, label = 'Prediction Interval')
for i in range(10):
y_tilde = gamma(a=shape, scale=np.exp(0.8 + X_pred[:, 1]*true_beta)/shape).rvs()
plt.scatter(X_pred[:, 1], y_tilde, s = 1, color = 'k')
plt.scatter(X_pred[:, 1], y_tilde, s = 1, color = 'k', label = 'New Data')
plt.legend()
Matematyka tego, co się dzieje
Nasze dane $y$ są dystrybuowane zgodnie z
$$ y\vert X \sim \mbox{Gamma}(\phi, \mu(x)/\phi) $$
Przynajmniej myślę, że to poprawna parametryzacja Gammy, nigdy nie mogę tego zrobić dobrze. W każdym razie, zakładając, że do modelu używamy łącza dziennika, oznacza to
$$ \mu(x) = \exp(X\beta)$$
Rzecz w tym, że nigdy nie wiemy $\beta$, mamy tylko $\hat{\beta}$ponieważ musimy oszacować parametry modelu. Parametry są zatem zmienną losową (ponieważ różne dane mogą dawać różne parametry). Teoria mówi, że przy wystarczającej ilości danych możemy rozważyć
$$ \hat{\beta} \sim \mbox{Normal}(\beta, \Sigma) $$
a trochę więcej teorii mówi, że podłączając nasze oszacowanie dla $\beta$ i $\Sigma$powinno być wystarczająco dobre. Pozwolić$\tilde{y}\vert X$ być danymi, które mogę zobaczyć w przypadku obserwacji ze zmiennymi towarzyszącymi $X$. Gdybym mógł, naprawdę bym liczył
$$ \tilde{y} \vert X \sim \int p(y\vert X,\beta)p (\beta) \, d \beta $$
a następnie weź kwantyle tej dystrybucji. Ale ta całka jest naprawdę trudna, więc zamiast tego po prostu ją przybliżamy, symulując z$p(\beta)$ (rozkład normalny) i przekazując wszystko, do czego symulowaliśmy $p(y\vert X, \beta)$ (w tym przypadku rozkład gamma).
Teraz zdaję sobie sprawę, że byłem tutaj dość szybki i luźny, więc jeśli któryś z czytelników chce włożyć trochę więcej rygoru w moje wyjaśnienie, daj mi znać w komentarzu, a ja to posprzątam. Myślę, że to powinno wystarczyć, aby pokazać OP, jak to działa.