StatsModels的置信度和预测间隔

Question

我这样做linear regression有StatsModels:

import numpy as np
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std

n = 100

x = np.linspace(0, 10, n)
e = np.random.normal(size=n)
y = 1 + 0.5*x + 2*e
X = sm.add_constant(x)

re = sm.OLS(y, X).fit()
print(re.summary())

prstd, iv_l, iv_u = wls_prediction_std(re)

我的问题是,iv_l和iv_u为上,下置信区间或预测区间？

我如何得到别人？

我需要所有点的置信度和预测间隔,做一个情节.

Answer 1

更新见第二个答案,这是最近的.一些模型和结果类现在具有get_prediction提供附加信息的方法,包括预测平均值的预测间隔和/或置信区间.

老答案:

iv_l并iv_u给出每个点的预测间隔的限制.

预测间隔是观察的置信区间,包括误差的估计.

我认为,平均预测的置信区间尚不可用statsmodels.(实际上,拟合值的置信区间隐藏在influence_outlier的summary_table中,但我需要验证这一点.)

适用于statsmodel的预测方法在TODO列表中.

加成

OLS存在置信区间,但访问有点笨拙.

要在运行脚本后包含:

from statsmodels.stats.outliers_influence import summary_table

st, data, ss2 = summary_table(re, alpha=0.05)

fittedvalues = data[:, 2]
predict_mean_se  = data[:, 3]
predict_mean_ci_low, predict_mean_ci_upp = data[:, 4:6].T
predict_ci_low, predict_ci_upp = data[:, 6:8].T

# Check we got the right things
print np.max(np.abs(re.fittedvalues - fittedvalues))
print np.max(np.abs(iv_l - predict_ci_low))
print np.max(np.abs(iv_u - predict_ci_upp))

plt.plot(x, y, 'o')
plt.plot(x, fittedvalues, '-', lw=2)
plt.plot(x, predict_ci_low, 'r--', lw=2)
plt.plot(x, predict_ci_upp, 'r--', lw=2)
plt.plot(x, predict_mean_ci_low, 'r--', lw=2)
plt.plot(x, predict_mean_ci_upp, 'r--', lw=2)
plt.show()

这应该给出与SAS相同的结果,http://jpktd.blogspot.ca/2012/01/nice-thing-about-seeing-zeros.html

Answer 2

对于测试数据,您可以尝试使用以下内容.

predictions = result.get_prediction(out_of_sample_df)
predictions.summary_frame(alpha=0.05)

我发现summary_frame()方法埋在这里,你可以找到get_prediction()方法在这里.您可以通过修改"alpha"参数来更改置信区间和预测区间的显着性级别.

我在这里发布这个是因为这是在寻找信心和预测间隔的解决方案时出现的第一篇文章 - 尽管这与测试数据相关.

以下是使用此方法获取模型,新数据和任意分位数的函数:

def ols_quantile(m, X, q):
  # m: OLS model.
  # X: X matrix.
  # q: Quantile.
  #
  # Set alpha based on q.
  a = q * 2
  if q > 0.5:
    a = 2 * (1 - q)
  predictions = m.get_prediction(X)
  frame = predictions.summary_frame(alpha=a)
  if q > 0.5:
    return frame.obs_ci_upper
  return frame.obs_ci_lower

Answer 3

通过时间序列结果，您可以使用该方法获得更加平滑的绘图get_forecast()。时间序列的示例如下：

# Seasonal Arima Modeling, no exogenous variable
model = SARIMAX(train['MI'], order=(1,1,1), seasonal_order=(1,1,0,12), enforce_invertibility=True)

results = model.fit()

results.summary()

下一步是进行预测，这会生成置信区间。

# make the predictions for 11 steps ahead
predictions_int = results.get_forecast(steps=11)
predictions_int.predicted_mean

这些可以放入数据框中，但需要一些清理：

# get a better view
predictions_int.conf_int()

连接数据框，但清理标题

conf_df = pd.concat([test['MI'],predictions_int.predicted_mean, predictions_int.conf_int()], axis = 1)

conf_df.head()

然后我们重命名这些列。

conf_df = conf_df.rename(columns={0: 'Predictions', 'lower MI': 'Lower CI', 'upper MI': 'Upper CI'})
conf_df.head()

制定情节。

# make a plot of model fit
# color = 'skyblue'

fig = plt.figure(figsize = (16,8))
ax1 = fig.add_subplot(111)


x = conf_df.index.values


upper = conf_df['Upper CI']
lower = conf_df['Lower CI']

conf_df['MI'].plot(color = 'blue', label = 'Actual')
conf_df['Predictions'].plot(color = 'orange',label = 'Predicted' )
upper.plot(color = 'grey', label = 'Upper CI')
lower.plot(color = 'grey', label = 'Lower CI')

# plot the legend for the first plot
plt.legend(loc = 'lower left', fontsize = 12)


# fill between the conf intervals
plt.fill_between(x, lower, upper, color='grey', alpha='0.2')

plt.ylim(1000,3500)

plt.show()