python: 模型的统计信息
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
from statsmodels.graphics.mosaicplot import mosaic
%pylab inline
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号 '-' 显示为方块的问题
载入和分析数据¶
def readData(path):
"""
使用pandas读取数据
"""
data = pd.read_csv(path)
cols = ["age", "education_num", "capital_gain", "capital_loss", "hours_per_week", "label"]
return data[cols]
data = readData('dataset/adult.data')
data.hist(rwidth=0.9, grid=False, figsize=(8, 8), alpha=0.6) # 各个变量的直方图
plt.show()
data.head()
data['label_code'] = pd.Categorical(data.label).codes # 将文字变量转化为数字变量
data[['label', 'label_code']].head(8)
def visualData(df):
"""
画直方图，直观了解数据
"""
df.hist(
rwidth=0.9, grid=False, figsize=(8, 8), alpha=0.6, color="grey")
plt.show(block=False)
df = data[["age", "education_num", "capital_gain", "capital_loss", "hours_per_week", "label_code"]]
visualData(df)
data.describe() # 仅仅显示数值型变量的基本统计信息，若要显示全部变量的统计信息只需传入参数 `include='all'`
交叉报表¶
交叉报表用来描述两个变量之间的关系。
- 按四分位数划分为 $4$ 个区间
qeducation_num = pd.qcut(data['education_num'], [.25, .5, .75, 1])
qeducation_num
- 计算
education_num和label的交叉报表
cross1 = pd.crosstab(data['label'], qeducation_num)
cross1
将 cross1 图表化
props = lambda key: {"color": "0.45"} if ' >50K' in key else {"color": "#C6E2FF"}
mosaic(cross1.stack(), gap=0.05,properties=props, title='label 和 education_num 的交叉报表')
plt.show()
def trainModel(data):
"""
搭建逻辑回归模型，并训练模型
"""
formula = "label_code ~ age + education_num + capital_gain + capital_loss + hours_per_week"
model = sm.Logit.from_formula(formula, data=data)
re = model.fit()
return re
def modelSummary(re):
"""
分析逻辑回归模型的统计性质
"""
# 整体统计分析结果
print(re.summary())
# 用f test检验education_num的系数是否显著
print("检验假设education_num的系数等于0：")
print(re.f_test("education_num=0"))
# 用f test检验两个假设是否同时成立
print("检验假设education_num的系数等于0.32和hours_per_week的系数等于0.04同时成立：")
print(re.f_test("education_num=0.32, hours_per_week=0.04"))
def interpretModel(re):
"""
理解模型结果
参数
----
re ：BinaryResults，训练好的逻辑回归模型
"""
conf = re.conf_int()
conf['OR'] = re.params
# 计算各个变量对事件发生比的影响
# conf里面的三列，分别对应着估计值的下界、上界和估计值本身
conf.columns = ['2.5%', '97.5%', 'OR']
print("各个变量对事件发生比的影响：")
print(np.exp(conf))
# 计算各个变量的边际效应
print("各个变量的边际效应：")
print(re.get_margeff(at="overall").summary())
def makePrediction(re, testSet, alpha=0.5):
"""
使用训练好的模型对测试数据做预测
"""
# 关闭pandas有关chain_assignment的警告
pd.options.mode.chained_assignment = None
# 计算事件发生的概率
testSet["prob"] = re.predict(testSet)
print("事件发生概率（预测概率）大于0.6的数据个数：")
print(testSet[testSet["prob"] > 0.6].shape[0]) # 输出值为576
print("事件发生概率（预测概率）大于0.5的数据个数：")
print(testSet[testSet["prob"] > 0.5].shape[0]) # 输出值为834
# 根据预测的概率，得出最终的预测
testSet["pred"] = testSet.apply(lambda x: 1 if x["prob"] > alpha else 0, axis=1)
return testSet
def evaluation(re):
"""
计算预测结果的查准查全率以及f1
参数
----
re ：DataFrame，预测结果，里面包含两列：真实值‘lable_code’、预测值‘pred’
"""
bins = np.array([0, 0.5, 1])
label = re["label_code"]
pred = re["pred"]
tp, fp, fn, tn = np.histogram2d(label, pred, bins=bins)[0].flatten()
precision = tp / (tp + fp) # 0.951
recall = tp / (tp + fn) # 0.826
f1 = 2 precision recall / (precision + recall) # 0.884
print("查准率: %.3f, 查全率: %.3f, f1: %.3f" % (precision, recall, f1))
# 将数据分为训练集和测试集
trainSet, testSet = train_test_split(data, test_size=0.2, random_state=2310)
# 训练模型并分析模型效果
re = trainModel(trainSet)
modelSummary(re)
interpretModel(re)
re = makePrediction(re, testSet)
evaluation(re)
