Bagaimana cara menghitung interval prediksi di GLM (Gamma) / TweedieRegression dengan Python?
Saya telah memeriksa banyak sumber dari web tentang melakukan interval prediksi, terutama dalam fungsi GLM. Salah satu pendekatannya adalah tentang Prediksi Interval untuk Machine Learninghttps://machinelearningmastery.com/prediction-intervals-for-machine-learning/dari Jason Brownlee. Namun, metodenya menargetkan regresi linier, dan mungkin tidak sesuai dengan GLM (Gamma) sampai derajat tertentu. Pendekatan lain yang saya temukan adalah menggunakan metode bootstrap untuk melakukan interval prediksi. Namun, perhitungannya sangat memakan waktu, dan memori komputer saya mati ketika menjalankan fungsi dari artikel tersebuthttps://saattrupdan.github.io/2020-03-01-bootstrap-prediction/. Saya bingung bagaimana melakukan prediksi interval dengan cara yang sesuai di GLM (kemungkinan besar Gamma) dengan Python, bukan di R. Saya telah menemukan paket terkait di R, tetapi saya tidak ingin menggunakan R untuk melakukan interval. Informasi terkait lainnya yang saya temukan dari web adalah Gamma GLM - Turunkan interval prediksi untuk yang baru x_i: Gamma GLM - Turunkan interval prediksi untuk x_i baru .
Jawaban
Agak terlibat, tetapi itu harus bisa dilakukan.
Seperti yang dikatakan posting itu, untuk mendapatkan interval prediksi, Anda harus mengintegrasikan ketidakpastian dalam koefisien. Itu sulit untuk dilakukan secara analitis, tetapi kita dapat mensimulasikannya. Berikut beberapa data regresi 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()
Sekarang, saya akan menyesuaikan regresi gamma ke data ini
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
Selama saya memiliki cukup data, kita dapat membuat perkiraan normal untuk distribusi sampling koefisien.
Mean dan kovarian dapat diperoleh dari ringkasan model.
beta_samp_mean = gamma_results.params
beta_samp_cov = gamma_results.cov_params()
dispersion = gamma_results.scale
Sekarang, ini hanya masalah pengambilan sampel data palsu menggunakan perkiraan ini dan mengambil jumlah.
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)
Mudah untuk memplot interval prediksi
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()
Matematika tentang apa yang sedang terjadi
Data kami $y$ didistribusikan menurut
$$ y\vert X \sim \mbox{Gamma}(\phi, \mu(x)/\phi) $$
Setidaknya saya pikir itu adalah parameterisasi yang benar dari Gamma, saya tidak pernah bisa melakukannya dengan benar. Bagaimanapun, dengan asumsi kami menggunakan tautan log untuk model, ini berarti
$$ \mu(x) = \exp(X\beta)$$
Masalahnya, kita tidak pernah tahu $\beta$, kami hanya mendapatkan $\hat{\beta}$karena kita harus memperkirakan parameter model. Parameter dengan demikian merupakan variabel acak (karena data yang berbeda dapat menghasilkan parameter yang berbeda). Teori mengatakan bahwa dengan data yang cukup, kita bisa mempertimbangkan
$$ \hat{\beta} \sim \mbox{Normal}(\beta, \Sigma) $$
dan beberapa teori lagi mengatakan bahwa memasukkan perkiraan kami untuk $\beta$ dan $\Sigma$seharusnya cukup baik. Membiarkan$\tilde{y}\vert X$ menjadi data yang mungkin saya lihat untuk observasi dengan kovariat $X$. Jika saya bisa, saya akan benar-benar menghitung
$$ \tilde{y} \vert X \sim \int p(y\vert X,\beta)p (\beta) \, d \beta $$
dan kemudian ambil kuantil dari distribusi ini. Tapi integral ini sangat sulit, jadi kita hanya memperkirakannya dengan mensimulasikan dari$p(\beta)$ (distribusi normal) dan meneruskan apa pun yang kami simulasi $p(y\vert X, \beta)$ (dalam hal ini, distribusi gamma).
Sekarang, saya menyadari bahwa saya telah cukup cepat dan longgar di sini, jadi jika ada pembaca yang ingin memberikan penjelasan saya sedikit lebih keras, beri tahu saya di komentar dan saya akan membersihkannya. Saya pikir ini seharusnya cukup baik untuk memberi OP gambaran tentang cara kerjanya.