PythonでGLM（ガンマ）/ TweedieRegressionで予測区間を計算する方法は？
特にGLM関数で、予測間隔の実行についてWebから多くのソースを確認しました。アプローチの1つは、機械学習の予測区間に関するものですhttps://machinelearningmastery.com/prediction-intervals-for-machine-learning/ジェイソンブラウンリーから。ただし、彼の方法は線形回帰を対象としているため、GLM（ガンマ）にはある程度適切ではない可能性があります。私が見つけたもう1つのアプローチは、ブートストラップ法を使用して予測区間を実行することです。ただし、計算には非常に時間がかかり、記事の関数を実行するとコンピューターのメモリが停止しました。https://saattrupdan.github.io/2020-03-01-bootstrap-prediction/。RではなくPythonのGLM（おそらくガンマ）で適切な方法で予測区間を実行する方法がわかりません。Rで関連パッケージを見つけましたが、Rを使用して区間を実行したくありません。Webから見つけたもう1つの関連情報は、ガンマGLM-新しいx_iの予測区間の導出x_i：ガンマGLM-新しいx_iの予測区間の導出です。
回答
少し複雑ですが、実行可能である必要があります。
その投稿が言うように、予測区間を取得するには、係数の不確実性を積分する必要があります。これを分析的に行うのは難しいですが、代わりにシミュレートすることができます。これがいくつかのガンマ回帰データです
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)$ （この場合、ガンマ分布）。
さて、私はここでかなり速くて緩んでいることに気づきました。読者が私の説明をもう少し厳密にしたい場合は、コメントで知らせてください。クリーンアップします。これは、これがどのように機能するかをOPに理解させるのに十分なはずだと思います。