statsmodels︱python常规统计模型库---素质

之前看sklearn线性模型没有R方，F检验，回归系数T检验等指标，于是看到了statsmodels这个库，看着该库输出的结果真是够怀念的。。

文章目录1 安装2 相关模型介绍2.1 线性模型2.2 离散选择模型(Discrete Choice Model, DCM)2.3 非参数统计2.4 广义线性模型 - Generalized Linear Models2.5 稳健回归——Robust Regression2.6 广义估计方程2.7 方差分析2.8 时间序列分析——Time Series Analysis2.9 空间计量必备：状态空间模型——State space models2.10 多元统计模型——因子/主成分分析3 相关模型demo3.1 线性回归模型3.2 广义线性模型——GLM3.3 稳健回归4 其他4.1 模型结果如何CSV导出？4.2 画模型图以及保存4.3 快速获取模型输出参数：P检验、F检验、P统计量1 安装pip install statsmodels

不过有可能会报错：

ImportError: cannot import name 'factorial' from 'scipy.misc' (E:Anaconda3.7libsite-packagesscipymisc__init__.py)

是跟scipy版本不匹配，笔者是删掉之前的pip uninstall statsmodels，再重新安装了一下就好了：

pip install --pre statsmodels -i https://pypi.tuna.tsinghua.edu.cn/simple2 相关模型介绍

相关文档可见：https://www.statsmodels.org/stable/examples/index.html

包含的模型有：

2.1 线性模型2.2 离散选择模型(Discrete Choice Model, DCM)

参考：离散选择模型(Discrete Choice Model, DCM)简介——之一

离散选择模型（Discrete Choice Model, DCM）在经济学领域和社会学领域都有广泛的应用。 例如，消费者在购买汽车的时候通常会比较几个不同的品牌，如福特、本田、大众，等等。 如果将消费者选择福特汽车记为Y=1，选择本田汽车记为Y=2，选择大众汽车记为Y=3；那么在研究消费者选择何种汽车品牌的时候，由于因变量不是一个连续的变量（Y=1, 2, 3），传统的线性回归模型就有一定的局限（见DCM系列文章第2篇）。 再比如，在交通安全研究领域，通常将交通事故的严重程度划分为3大类：

（1）仅财产损失（Property Damage Only, PDO），（2）受伤（Injury），（3）死亡（Fatality）； 在研究各类因素（如道路坡度、弯道曲率等、车龄、光照、天气条件等）对事故严重程度的影响的时候，由于因变量（事故严重程度）是一个离散变量（仅3个选项），使用离散选择模型可以提供一个有效的建模途径。 2.3 非参数统计2.4 广义线性模型 - Generalized Linear Models2.5 稳健回归——Robust Regression2.6 广义估计方程2.7 方差分析2.8 时间序列分析——Time Series Analysis2.9 空间计量必备：状态空间模型——State space models2.10 多元统计模型——因子/主成分分析3 相关模型demo3.1 线性回归模型

可参考：https://www.statsmodels.org/stable/examples/notebooks/generated/ols.html

# 线性模型 import statsmodels.api as sm import numpy as np x = np.linspace(0,10,100) y = 3*x + np.random.randn()+ 10 # Fit and summarize OLS model X = sm.add_constant(x) mod = sm.OLS(y,X) result = mod.fit() print('Parameters: ', result .params) print('Standard errors: ', result .bse) print('Predicted values: ', result .predict()) print(result.summary()) # 预测数据 print(result.predict(X[:5]))

输出结果超级熟悉。

result.params是回归系数result.summary()把模型相关系数都打印出来 其中，预测的时候，如果不给入参数result.predict()，则默认是X3.2 广义线性模型——GLM

参考：https://www.statsmodels.org/stable/examples/notebooks/generated/glm.html

import statsmodels.formula.api as smf star98 = sm.datasets.star98.load_pandas().data formula = 'SUCCESS ~ LOWINC + PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF' dta = star98[['NABOVE', 'NBELOW', 'LOWINC', 'PERASIAN', 'PERBLACK', 'PERHISP', 'PCTCHRT', 'PCTYRRND', 'PERMINTE', 'AVYRSEXP', 'AVSALK', 'PERSPENK', 'PTRATIO', 'PCTAF']].copy() endog = dta['NABOVE'] / (dta['NABOVE'] + dta.pop('NBELOW')) del dta['NABOVE'] dta['SUCCESS'] = endog mod1 = smf.glm(formula=formula, data=dta, family=sm.families.Binomial()).fit() mod1.summary() mod1.predict(dta)

formula是常规的公式，其中所有X/Y数据都放在一个dataframe之中。

print('Total number of trials:', data.endog[0].sum()) print('Parameters: ', res.params) print('T-values: ', res.tvalues)

包括了回归系数，T检验值

3.3 稳健回归

参考：https://www.statsmodels.org/stable/examples/notebooks/generated/robust_models_0.html

nsample = 50 x1 = np.linspace(0, 20, nsample) X = np.column_stack((x1, (x1-5)**2)) X = sm.add_constant(X) sig = 0.3 # smaller error variance makes OLS<->RLM contrast bigger beta = [5, 0.5, -0.0] y_true2 = np.dot(X, beta) y2 = y_true2 + sig*1. * np.random.normal(size=nsample) y2[[39,41,43,45,48]] -= 5 # add some outliers (10% of nsample) X2 = X[:,[0,1]] res2 = sm.OLS(y2, X2).fit() print(res2.params) print(res2.bse) resrlm2 = sm.RLM(y2, X2).fit() print(resrlm2.params) print(resrlm2.bse) print(resrlm2.summary())4 其他4.1 模型结果如何CSV导出？

可以通过as_csv()将模型导出

resrlm2 = sm.RLM(y, x).fit() resrlm2.summary() with open( 'model_rlm.csv', 'w') as fh: fh.write(resrlm2.summary().as_csv())

不过导出的格式比较奇怪：

4.2 画模型图以及保存import statsmodels.api as sm import numpy as np import matplotlib.pyplot as plt # 准备数据 x = np.linspace(0,10,100) y = 3*x + np.random.randn()+ 10 # Fit and summarize OLS model res = sm.OLS(y,x).fit() print(res.params) print(res.summary()) # 稳健回归 resrlm = sm.RLM(y, x).fit() # 画图 fig, ax = plt.subplots(figsize=(8,6)) ax.plot(x, y, 'o', label="truey ") ax.plot(x, res.predict(), 'o', label="ols") # res2.predict(X2) == res2.predict() ax.plot(x, resrlm.predict(), 'b-', label="rlm")# resrlm2.predict(X2) == resrlm2.predict() legend = ax.legend(loc="best") # 图保存 plt.savefig( 'image.jpg')4.3 快速获取模型输出参数：P检验、F检验、P统计量def get_model_param(res2,name = 'all'): model_param_dict = {'name':name, # 模型的名字 'rsquared':res2.rsquared, # R方 'fvalue':res2.fvalue, # F值，整个模型 'f_pvalue':res2.f_pvalue, # P值，整个模型 'params':res2.params[0], # 回归系数 'pvalues':res2.pvalues[0], # 回归系数 P检验 0.000 'tvalues':res2.tvalues[0]} # 回归系数 T检验 276.571 return model_param_dict ---来自腾讯云社区的---素质