你可以使用
增量法来找到预测概率的近似方差。也就是说,
var(proba) = np.dot(np.dot(gradient.T, cov), gradient)
其中gradient
是模型系数预测概率的导数向量,cov
是系数协方差矩阵。
Delta方法已被证明在所有最大似然估计中渐近地有效。但是,如果您的训练样本很小,则渐近方法可能效果不佳,您应该考虑引入自助法。
以下是将Delta方法应用于逻辑回归的示例:
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
np.random.seed(1)
x = np.arange(100)
y = (x * 0.5 + np.random.normal(size=100,scale=10)>30)
X = sm.add_constant(x)
model = sm.Logit(y, X).fit()
proba = model.predict(X)
cov = model.cov_params()
gradient = (proba * (1 - proba) * X.T).T
std_errors = np.array([np.sqrt(np.dot(np.dot(g, cov), g)) for g in gradient])
c = 1.96
upper = np.maximum(0, np.minimum(1, proba + std_errors * c))
lower = np.maximum(0, np.minimum(1, proba - std_errors * c))
plt.plot(x, proba)
plt.plot(x, lower, color='g')
plt.plot(x, upper, color='g')
plt.show()
它绘制了以下漂亮的图片:![enter image description here](https://istack.dev59.com/7oZ4D.webp)
对于您的示例,代码将是:
proba = logit.predict(age_range_poly)
cov = logit.cov_params()
gradient = (proba * (1 - proba) * age_range_poly.T).T
std_errors = np.array([np.sqrt(np.dot(np.dot(g, cov), g)) for g in gradient])
c = 1.96
upper = np.maximum(0, np.minimum(1, proba + std_errors * c))
lower = np.maximum(0, np.minimum(1, proba - std_errors * c))
plt.plot(age_range_poly[:, 1], proba)
plt.plot(age_range_poly[:, 1], lower, color='g')
plt.plot(age_range_poly[:, 1], upper, color='g')
plt.show()
然后它将呈现以下图片
![enter image description here](https://istack.dev59.com/ld7Qs.webp)
看起来很像一条囊括了一只大象的巨蟒。
您可以将其与引导估计进行比较:
preds = []
for i in range(1000):
boot_idx = np.random.choice(len(age), replace=True, size=len(age))
model = sm.Logit(wage['wage250'].iloc[boot_idx], age[boot_idx]).fit(disp=0)
preds.append(model.predict(age_range_poly))
p = np.array(preds)
plt.plot(age_range_poly[:, 1], np.percentile(p, 97.5, axis=0))
plt.plot(age_range_poly[:, 1], np.percentile(p, 2.5, axis=0))
plt.show()
![enter image description here](https://istack.dev59.com/GQZT3.webp)
差分法和 bootstrap 的结果看起来非常相似。
然而,这本书的作者采用了第三种方法。他们利用了以下事实:
proba = np.exp(np.dot(x, params)) / (1 + np.exp(np.dot(x, params)))
并计算线性部分的置信区间,然后再使用 logit 函数进行转换。
xb = np.dot(age_range_poly, logit.params)
std_errors = np.array([np.sqrt(np.dot(np.dot(g, cov), g)) for g in age_range_poly])
upper_xb = xb + c * std_errors
lower_xb = xb - c * std_errors
upper = np.exp(upper_xb) / (1 + np.exp(upper_xb))
lower = np.exp(lower_xb) / (1 + np.exp(lower_xb))
plt.plot(age_range_poly[:, 1], upper)
plt.plot(age_range_poly[:, 1], lower)
plt.show()
所以他们得到了分歧区间:
![enter image description here](https://istack.dev59.com/haApj.webp)
这些方法产生如此不同的结果,是因为它们假定分布正常的“不同事物”(预测的概率和对数几率)。也就是说,Delta方法假定预测的概率是正常的,在书中,对数几率是正常的。实际上,它们在有限样本中都不是正态的,它们都在无限样本中收敛于正态,但它们的方差同时收敛于零。最大似然估计对重新参数化不敏感,但它们的估计分布是敏感的,这就是问题所在。