最大似然估计（通用模型）¶
本教程解释了如何在statsmodels中快速实现新的最大似然模型。我们给出两个例子：
二元因变量的概率模型
计数数据的负二项模型
类通过提供诸如自动数值微分和统一接口等工具来简化过程，这些接口可以与优化函数一起使用。使用，用户只需“插入”对数似然函数即可轻松拟合新的MLE模型。
示例 1: Probit 模型¶
[1]:
import numpy as np
from scipy import stats
import statsmodels.api as sm
from statsmodels.base.model import GenericLikelihoodModel
The Spector 数据集随 statsmodels 一起分发。您可以像这样访问因变量 (endog) 的值向量和回归变量矩阵 (exog)：
[2]:
data = sm.datasets.spector.load_pandas()
exog = data.exog
endog = data.endog
print(sm.datasets.spector.NOTE)
print(data.exog.head())
::
Number of Observations - 32
Number of Variables - 4
Variable name definitions::
Grade - binary variable indicating whether or not a student's grade
improved. 1 indicates an improvement.
TUCE - Test score on economics test
PSI - participation in program
GPA - Student's grade point average
GPA TUCE PSI
0 2.66 20.0 0.0
1 2.89 22.0 0.0
2 3.28 24.0 0.0
3 2.92 12.0 0.0
4 4.00 21.0 0.0
然后，我们向回归矩阵添加一个常数：
[3]:
exog = sm.add_constant(exog, prepend=True)
要创建您自己的似然模型，您只需覆盖 loglike 方法。
[4]:
class MyProbit(GenericLikelihoodModel):
def loglike(self, params):
exog = self.exog
endog = self.endog
q = 2 * endog - 1
return stats.norm.logcdf(q*np.dot(exog, params)).sum()
估计模型并打印摘要：
[5]:
sm_probit_manual = MyProbit(endog, exog).fit()
print(sm_probit_manual.summary())
Optimization terminated successfully.
Current function value: 0.400588
Iterations: 292
Function evaluations: 494
MyProbit Results
==============================================================================
Dep. Variable: GRADE Log-Likelihood: -12.819
Model: MyProbit AIC: 33.64
Method: Maximum Likelihood BIC: 39.50
Date: Wed, 16 Oct 2024
Time: 18:27:05
No. Observations: 32
Df Residuals: 28
Df Model: 3
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -7.4523 2.542 -2.931 0.003 -12.435 -2.469
GPA 1.6258 0.694 2.343 0.019 0.266 2.986
TUCE 0.0517 0.084 0.617 0.537 -0.113 0.216
PSI 1.4263 0.595 2.397 0.017 0.260 2.593
==============================================================================
将您的Probit实现与statsmodels的“预制”实现进行比较：
[6]:
sm_probit_canned = sm.Probit(endog, exog).fit()
Optimization terminated successfully.
Current function value: 0.400588
Iterations 6
[7]:
print(sm_probit_canned.params)
print(sm_probit_manual.params)
const -7.452320
GPA 1.625810
TUCE 0.051729
PSI 1.426332
dtype: float64
[-7.45233176 1.62580888 0.05172971 1.42631954]
[8]:
print(sm_probit_canned.cov_params())
print(sm_probit_manual.cov_params())
const GPA TUCE PSI
const 6.464166 -1.169668 -0.101173 -0.594792
GPA -1.169668 0.481473 -0.018914 0.105439
TUCE -0.101173 -0.018914 0.007038 0.002472
PSI -0.594792 0.105439 0.002472 0.354070
[[ 6.46416763e+00 -1.16966612e+00 -1.01173185e-01 -5.94788994e-01]
[-1.16966612e+00 4.81472098e-01 -1.89134583e-02 1.05438216e-01]
[-1.01173185e-01 -1.89134583e-02 7.03758407e-03 2.47189357e-03]
[-5.94788994e-01 1.05438216e-01 2.47189357e-03 3.54069511e-01]]
请注意，GenericMaximumLikelihood 类提供了自动微分功能，因此我们不需要提供 Hessian 或 Score 函数来计算协方差估计。
示例 2: 计数数据的负二项回归¶
考虑一个用于计数数据的负二项回归模型，其对数似然（类型NB-2）函数表示为：
\[\mathcal{L}(\beta_j; y, \alpha) = \sum_{i=1}^n y_i ln
\left ( \frac{\alpha exp(X_i'\beta)}{1+\alpha exp(X_i'\beta)} \right ) -
\frac{1}{\alpha} ln(1+\alpha exp(X_i'\beta)) + ln \Gamma (y_i + 1/\alpha) - ln \Gamma (y_i+1) - ln \Gamma (1/\alpha)\]
使用回归变量矩阵 \(X\)，系数向量 \(\beta\)，以及负二项分布异质性参数 \(\alpha\)。
使用来自 scipy 的 nbinom 分布，我们可以简单地将这个似然写为：
[9]:
import numpy as np
from scipy.stats import nbinom
[10]:
def _ll_nb2(y, X, beta, alph):
mu = np.exp(np.dot(X, beta))
size = 1/alph
prob = size/(size+mu)
ll = nbinom.logpmf(y, size, prob)
return ll
新模型类¶
我们创建一个继承自 GenericLikelihoodModel 的新模型类：
[11]:
from statsmodels.base.model import GenericLikelihoodModel
[12]:
class NBin(GenericLikelihoodModel):
def __init__(self, endog, exog, **kwds):
super(NBin, self).__init__(endog, exog, **kwds)
def nloglikeobs(self, params):
alph = params[-1]
beta = params[:-1]
ll = _ll_nb2(self.endog, self.exog, beta, alph)
return -ll
def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):
# we have one additional parameter and we need to add it for summary
self.exog_names.append('alpha')
if start_params == None:
# Reasonable starting values
start_params = np.append(np.zeros(self.exog.shape[1]), .5)
# intercept
start_params[-2] = np.log(self.endog.mean())
return super(NBin, self).fit(start_params=start_params,
maxiter=maxiter, maxfun=maxfun,
**kwds)
需要注意的两点：
nloglikeobs: 该函数应返回数据集中每个观测值的负对数似然函数的一次评估（即 endog/X 矩阵的行）。start_params: 需要提供一个一维的初始值数组。该数组的大小决定了优化中将使用的参数数量。
就是这样！你完成了！
使用示例¶
数据集 Medpar 以CSV格式托管在 Rdatasets 仓库。我们使用 Pandas 库 中的 read_csv 函数将数据加载到内存中。然后我们打印前几列：
[13]:
import statsmodels.api as sm
[14]:
medpar = sm.datasets.get_rdataset("medpar", "COUNT", cache=True).data
medpar.head()
[14]:
| los | hmo | white | died | age80 | type | type1 | type2 | type3 | provnum | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 4 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 30001 |
| 1 | 9 | 1 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 30001 |
| 2 | 3 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 30001 |
| 3 | 9 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 30001 |
| 4 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 30001 |
我们感兴趣的模型有一个非负整数向量作为因变量（los），以及5个回归变量：Intercept、type2、type3、hmo、white。
对于估计，我们需要创建两个变量来保存我们的回归变量和结果变量。这些可以是ndarrays或pandas对象。
[15]:
y = medpar.los
X = medpar[["type2", "type3", "hmo", "white"]].copy()
X["constant"] = 1
然后，我们拟合模型并提取一些信息：
[16]:
mod = NBin(y, X)
res = mod.fit()
Optimization terminated successfully.
Current function value: 3.209014
Iterations: 805
Function evaluations: 1238
/Users/cw/baidu/code/fin_tool/github/statsmodels/venv/lib/python3.11/site-packages/statsmodels/base/model.py:2748: UserWarning: df_model + k_constant + k_extra differs from k_params
warnings.warn("df_model + k_constant + k_extra "
/Users/cw/baidu/code/fin_tool/github/statsmodels/venv/lib/python3.11/site-packages/statsmodels/base/model.py:2752: UserWarning: df_resid differs from nobs - k_params
warnings.warn("df_resid differs from nobs - k_params")
提取参数估计值、标准误差、p值、AIC等：
[17]:
print('Parameters: ', res.params)
print('Standard errors: ', res.bse)
print('P-values: ', res.pvalues)
print('AIC: ', res.aic)
Parameters: [ 0.2212642 0.70613942 -0.06798155 -0.12903932 2.31026565 0.44575147]
Standard errors: [0.0505926 0.07613046 0.05326097 0.06854122 0.06794674 0.01981542]
P-values: [1.22298424e-005 1.76975731e-020 2.01819150e-001 5.97474348e-002
2.14405153e-253 4.62692220e-112]
AIC: 9604.953205830161
像往常一样，您可以通过输入 dir(res) 来获取可用信息的完整列表。我们还可以查看估计结果的摘要。
[18]:
print(res.summary())
NBin Results
==============================================================================
Dep. Variable: los Log-Likelihood: -4797.5
Model: NBin AIC: 9605.
Method: Maximum Likelihood BIC: 9632.
Date: Wed, 16 Oct 2024
Time: 18:27:06
No. Observations: 1495
Df Residuals: 1490
Df Model: 4
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
type2 0.2213 0.051 4.373 0.000 0.122 0.320
type3 0.7061 0.076 9.275 0.000 0.557 0.855
hmo -0.0680 0.053 -1.276 0.202 -0.172 0.036
white -0.1290 0.069 -1.883 0.060 -0.263 0.005
constant 2.3103 0.068 34.001 0.000 2.177 2.443
alpha 0.4458 0.020 22.495 0.000 0.407 0.485
==============================================================================
测试¶
我们可以使用statsmodels实现的负二项模型来检查结果，该模型使用了解析得分函数和Hessian矩阵。
[19]:
res_nbin = sm.NegativeBinomial(y, X).fit(disp=0)
print(res_nbin.summary())
NegativeBinomial Regression Results
==============================================================================
Dep. Variable: los No. Observations: 1495
Model: NegativeBinomial Df Residuals: 1490
Method: MLE Df Model: 4
Date: Wed, 16 Oct 2024 Pseudo R-squ.: 0.01215
Time: 18:27:06 Log-Likelihood: -4797.5
converged: True LL-Null: -4856.5
Covariance Type: nonrobust LLR p-value: 1.404e-24
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
type2 0.2212 0.051 4.373 0.000 0.122 0.320
type3 0.7062 0.076 9.276 0.000 0.557 0.855
hmo -0.0680 0.053 -1.276 0.202 -0.172 0.036
white -0.1291 0.069 -1.883 0.060 -0.263 0.005
constant 2.3103 0.068 34.001 0.000 2.177 2.443
alpha 0.4457 0.020 22.495 0.000 0.407 0.485
==============================================================================
[20]:
print(res_nbin.params)
type2 0.221218
type3 0.706173
hmo -0.067987
white -0.129053
constant 2.310279
alpha 0.445748
dtype: float64
[21]:
print(res_nbin.bse)
type2 0.050592
type3 0.076131
hmo 0.053261
white 0.068541
constant 0.067947
alpha 0.019815
dtype: float64
或者我们可以将它们与使用R的MASS实现获得的结果进行比较：
url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/csv/COUNT/medpar.csv'
medpar = read.csv(url)
f = los~factor(type)+hmo+white
library(MASS)
mod = glm.nb(f, medpar)
coef(summary(mod))
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.31027893 0.06744676 34.253370 3.885556e-257
factor(type)2 0.22124898 0.05045746 4.384861 1.160597e-05
factor(type)3 0.70615882 0.07599849 9.291748 1.517751e-20
hmo -0.06795522 0.05321375 -1.277024 2.015939e-01
white -0.12906544 0.06836272 -1.887951 5.903257e-02
数值精度¶
The statsmodels 通用 MLE 和 R 参数估计在第四位小数之前是一致的。然而，标准误差仅在第二位小数之前一致。这种差异是由于我们对 Hessian 数值估计的不精确性造成的。在当前情况下，MASS 和 statsmodels 标准误差估计之间的差异在实质上是不相关的，但它突显了一个事实，即需要非常精确估计的用户可能并不总是希望在使用数值导数时依赖默认设置。在这种情况下，最好使用带有 LikelihoodModel 类的解析导数。
Last update:
Oct 16, 2024