पायथन में GLM (गामा) / TweedieRegression में भविष्यवाणी अंतराल की गणना कैसे करें?
मैंने विशेष रूप से जीएलएम फ़ंक्शन में प्रेडिटिटोन अंतराल के संचालन के बारे में जाले से बहुत स्रोत की जाँच की है। दृष्टिकोणों में से एक मशीन लर्निंग के लिए भविष्यवाणी अंतराल के बारे में हैhttps://machinelearningmastery.com/prediction-intervals-for-machine-learning/जेसन ब्राउनली से। हालाँकि, उनकी विधि रैखिक प्रतिगमन को लक्षित करती है, और यह कुछ डिग्री तक GLM (गामा) के लिए उपयुक्त नहीं हो सकता है। एक और दृष्टिकोण जो मैंने पाया है कि भविष्यवाणी अंतराल का संचालन करने के लिए बूटस्ट्रैपिंग विधि का उपयोग करना है। हालाँकि, गणना बहुत समय लेने वाली थी, और लेख से फ़ंक्शन चलाते समय मेरे कंप्यूटर की मेमोरी को मार दिया गया थाhttps://saattrupdan.github.io/2020-03-01-bootstrap-prediction/। मैं भ्रमित हूं कि आरएच के बजाय पायथन में जीएलएम (गामा सबसे संभवतया) में एक उपयुक्त तरीके से भविष्यवाणी अंतराल का संचालन कैसे करें। मुझे आर में संबंधित पैकेज मिला है, लेकिन मैं अंतराल का संचालन करने के लिए आर का उपयोग नहीं करना चाहता हूं। एक अन्य संबंधित जानकारी जो मुझे वेब से मिली वह है गामा GLM - नए के लिए व्युत्पन्न पूर्वानुमान अंतराल x_i: नए x_i के लिए गामा GLM - व्युत्पन्न पूर्वानुमान अंतराल ।
जवाब
इसकी एक बिट शामिल है, लेकिन यह योग्य होना चाहिए।
जैसा कि पोस्ट में कहा गया है, एक भविष्यवाणी अंतराल प्राप्त करने के लिए आपको गुणांकों में अनिश्चितता पर एकीकृत करना होगा। यह विश्लेषणात्मक रूप से करना कठिन है, लेकिन हम इसके बजाय अनुकरण कर सकते हैं। यहाँ कुछ गामा प्रतिगमन डेटा है
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()
अब, मैं इस डेटा के लिए गामा प्रतिगमन फिट करूँगा
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
जब तक मेरे पास पर्याप्त डेटा है, हम गुणांक के नमूना वितरण के लिए एक सामान्य सन्निकटन बना सकते हैं।
माध्य और सहसंयोजक मॉडल सारांश से प्राप्त किया जा सकता है।
beta_samp_mean = gamma_results.params
beta_samp_cov = gamma_results.cov_params()
dispersion = gamma_results.scale
अब, यह केवल इन अनुमानों का उपयोग करके और क्वांटिल लेने के लिए नकली डेटा का नमूना लेने की बात है।
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)
इसके बाद भविष्यवाणी अंतराल की साजिश करना आसान है
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()
जो चल रहा है उसका गणित
हमारा डेटा $y$ के अनुसार वितरित किए जाते हैं
$$ y\vert X \sim \mbox{Gamma}(\phi, \mu(x)/\phi) $$
कम से कम मुझे लगता है कि यह गामा का सही पैरामीटर है, मैं इसे कभी सही नहीं कर सकता। किसी भी मामले में, यह मानते हुए कि हम मॉडल के लिए लॉग लिंक का उपयोग करते हैं, इसका मतलब है
$$ \mu(x) = \exp(X\beta)$$
बात यह है, हम कभी नहीं जानते $\beta$, हम केवल प्राप्त करते हैं $\hat{\beta}$क्योंकि हमें मॉडल के मापदंडों का अनुमान लगाना होगा। इस प्रकार पैरामीटर एक यादृच्छिक चर हैं (क्योंकि विभिन्न डेटा विभिन्न मापदंडों को प्राप्त कर सकते हैं)। थ्योरी कहती है कि पर्याप्त डेटा के साथ, हम विचार कर सकते हैं
$$ \hat{\beta} \sim \mbox{Normal}(\beta, \Sigma) $$
और कुछ और सिद्धांत कहते हैं कि हमारे अनुमान में प्लगिंग $\beta$ तथा $\Sigma$काफी अच्छा होना चाहिए। लश्कर$\tilde{y}\vert X$ डेटा हो सकता है कि मैं कोविरेट के साथ टिप्पणियों के लिए देख सकता हूं $X$। अगर मैं कर सका, तो मैं वास्तव में गणना करूंगा
$$ \tilde{y} \vert X \sim \int p(y\vert X,\beta)p (\beta) \, d \beta $$
और फिर इस वितरण की मात्राएँ लें। लेकिन यह अभिन्नता वास्तव में कठिन है, इसलिए इसके बजाय हम केवल अनुकरण करके इसे अनुमानित करते हैं$p(\beta)$ (सामान्य वितरण) और जो कुछ भी हमने अनुकरण किया है उसे पारित करना $p(y\vert X, \beta)$ (इस मामले में, गामा वितरण)।
अब, मुझे एहसास हुआ कि मैं यहां काफी तेज और ढीला हूं, इसलिए यदि कोई भी पाठक मेरे स्पष्टीकरण में थोड़ी और सख्ती करना चाहते हैं, तो कृपया मुझे एक टिप्पणी में बताएं और मैं इसे साफ कर दूंगा। मुझे लगता है कि ओपी को यह अंदाजा लगाने के लिए यह काफी अच्छा होना चाहिए कि यह कैसे काम करता है।