我用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
是上下置信区间或预测区间?
我如何得到别人?
我需要所有点的置信度和预测间隔,来做一个情节。
更新见第二个答案,这是最近的。一些模型和结果类现在具有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相同的结果
对于测试数据,您可以尝试使用以下内容。
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
您可以使用我的repo(https://github.com/shahejokarian/regression-prediction-interval)中的Ipython笔记本中的LRPI()类来获取预测间隔。
您需要设置t值以获得预测值的所需置信区间,否则默认值为95%conf。间隔。
LRPI类使用sklearn.linear_model的LinearRegression,numpy和pandas库。
笔记本中也有一个例子。
当你需要单个分位数的精确结果时,summary_frame
和summary_table
运行良好,但不能很好地矢量化。这将提供预测区间的正常近似值(不是置信区间),并且适用于分位数矢量:
def ols_quantile(m, X, q):
# m: Statsmodels OLS model.
# X: X matrix of data to predict.
# q: Quantile.
#
from scipy.stats import norm
mean_pred = m.predict(X)
se = np.sqrt(m.scale)
return mean_pred + norm.ppf(q) * se
您可以根据statsmodel给出的结果和常态假设来计算它们。
以下是OLS和CI的平均值示例:
import statsmodels.api as sm
import numpy as np
from scipy import stats
#Significance level:
sl = 0.05
#Evaluate mean value at a required point x0. Here, at the point (0.0,2.0) for N_model=2:
x0 = np.asarray([1.0, 0.0, 2.0])# If you have no constant in your model, remove the first 1.0. For more dimensions, add the desired values.
#Get an OLS model based on output y and the prepared vector X (as in your notation):
model = sm.OLS(endog = y, exog = X )
results = model.fit()
#Get two-tailed t-values:
(t_minus, t_plus) = stats.t.interval(alpha = (1.0 - sl), df = len(results.resid) - len(x0) )
y_value_at_x0 = np.dot(results.params, x0)
lower_bound = y_value_at_x0 + t_minus*np.sqrt(results.mse_resid*( np.dot(np.dot(x0.T,results.normalized_cov_params),x0) ))
upper_bound = y_value_at_x0 + t_plus*np.sqrt(results.mse_resid*( np.dot(np.dot(x0.T,results.normalized_cov_params),x0) ))
你可以用输入结果,点x0和显着性水平sl包围一个很好的函数。
我现在不确定你是否可以将它用于WLS(),因为那里还有额外的事情发生。
参考:[D.C.蒙哥马利和E.A.啄。 “线性回归分析简介。”第4版。 Ed。,Wiley,1992]。