Как рассчитать интервал прогнозирования в GLM (Gamma) / TweedieRegression в Python?
Я проверил много источников в сети о проведении интервала предсказания, особенно в функции GLM. Один из подходов - интервалы прогнозирования для машинного обучения.https://machinelearningmastery.com/prediction-intervals-for-machine-learning/от Джейсона Браунли. Однако его метод нацелен на линейную регрессию, и в какой-то степени он может не подходить для GLM (гамма). Другой подход, который я нашел, - это использовать метод начальной загрузки для определения интервала прогнозирования. Однако вычисления были настолько трудоемкими, и память моего компьютера была убита при запуске функции из статьиhttps://saattrupdan.github.io/2020-03-01-bootstrap-prediction/. Я не понимаю, как правильно проводить интервал прогнозирования в GLM (скорее всего, гамма) в Python вместо R. Я нашел связанный пакет в R, но я не хочу использовать R для проведения интервала. Другая связанная информация, которую я нашел в Интернете, - это Gamma GLM - определение интервалов прогнозирования для новых x_i: Gamma 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 представление о том, как это работает.