stats模型的置信区间和预测区间

Question

我用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

对于测试数据，可以尝试使用以下方法。

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

我发现summary_frame（）方法隐藏在here中，您可以找到get_prediction（）方法here。通过修改“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 2

您可以使用我的repo（https://github.com/shahejokarian/regression-prediction-interval）中Ipython笔记本中的LRPI（）类来获取预测间隔。

您需要设置t值以获得预测值所需的置信区间，否则默认值为95%conf.interval。

LRPI类使用sklearn.linear_模型的LinearRegression、numpy和pandas库。

笔记本上也有一个例子。

Answer 3

更新请参阅第二个更新的答案。一些模型和结果类现在有一个get_prediction方法，该方法提供包括预测区间和/或预测平均值的置信区间在内的附加信息。

旧答案：

iv_l和iv_u给出了每个点的预测间隔的限制。

预测区间是观测值的置信区间，包括误差估计。

我认为，平均预测的置信区间在statsmodels中还不可用。 （实际上，拟合值的置信区间隐藏在影响异常值汇总表中，但我需要对此进行验证。）

正确的statsmodels预测方法在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