前言以及回归分析

统计建模

前言

这次想尝试一下统计建模，准备长三角建模的同时，加强一下自己的数据分析能力，学习的教材是张良均的《Python数据分析与挖掘实战》，这本书后面的实战练习相当不错，值得一做，书长下面这样，zlib里有pdf的，里面涉及的代码和数据在本书中均有给出获得方式。此学习笔记仅供参考，大部分内容和书中的内容相同，代码部分会有一些改进。

这本书第1~5章在讲基础部分，主要是数据的处理方法、python的使用以及一些常用的模型与方法，前四章内容我先前接触过，内容较为简单，这里就先略过，从第五章 挖掘建模开始学习。

第五章 挖掘建模

5.1 分类与预测

分类主要是预测分类标号（离散属性）， 而预测主要是建立连续值函数模型，预测给定自变量对应的因变量的值。

5.1.1 实现过程

5.1.2 常用算法

5.1.3 回归分析

回归分析：通过建立模型来研究变量之间相互关系的密切程度、结构状态及进行模型预测的有效工具。研究的范围大致如下：

在数据挖掘环境下，自变量与因变量具有相关关系，自变量的值是已知的，因变量是要预测的。

常用的回归模型见下表：

线性回归模型是相对简单的回归模型，但是通常因变量和自变量之间呈现出某种曲线关系，这就需要建立非线性回归模型。

1. Logistic 回归分析介绍

Logistic 回归属于概率型非线性回归，分为二分类和多分类的回归模型。

（ 1 ）Logistic 函数

Logistic 回归模型中的因变量只有 1-0(如“是”和“否”、“发生”和“不发生”)两种取值。假设在\(p\)个独立自变量\(x_1,x_2,\cdots,x_p\)的作用下，记\(y\)取 1 的概率是\(p=P(y=1\mid X)\) , 取 0 的概率是 1\(-p\) ,取 1 和取 0 的概率之比为\(\frac p{1-p}\),称为事件的优势比(odds),对odds 取自然对数即得 Logistic 变换 \(\mathrm{Logit}(p)=\ln\left(\frac p{1-p}\right)\circ\)

令 \(\mathrm{Logit}(p)=\ln\left(\frac p{1-p}\right)=z\), 则\(p=\frac1{1+\mathrm{e}^{-z}}\)即为 Logistic 函数，如图所示。

（ 2 ） Logistic 回归模型

Logistic 回归模型是建立 \(\ln\left(\frac p{1-p}\right)\)与自变量的线性回归模型。Logistic 回归模型为下式：

\[\ln\biggl(\frac{p}{1-p}\biggr)=\beta_0+\beta_1x_1+\cdots+\beta_px_p+\varepsilon \]

因为\(\ln\left(\frac p{1-n}\right)\)的取值范围是\((-\infty,+\infty)\),这样，自变量 \(x_1,x_2,\cdots,x_p\) 可在任意范围内取值。记\(g(x)=\beta_0+\beta_1x_1+\cdots+\beta_px_p\),得到式\((5-2)\)和式\((5-3)\)。

\[p=P(y=1\mid X)=\frac{1}{1+\mathrm{e}^{-g(x)}} \tag{5-2} \]

\[1-p=P(y=0\mid X)=1-\frac{1}{1+\mathrm{e}^{-g(x)}}=\frac{1}{1+\mathrm{e}^{g(x)}}\tag{5-3} \]

（ 3 ）Logistic 回归模型的解释

\[\frac p{1-p}=\mathrm{e}^{\beta_0+\beta_1x_1+\cdots+\beta_px_p+\varepsilon}\tag{5-4} \]

\(\beta_0:\)在没有自变量，即\(x_1,x_2,\cdots,x_p\)全部取 0 时，\(y=1\)与\(y=0\)发生概率之比的自然对数；

\(\beta_1:\)某自变量\(x_i\)变化时，即\(x_i=1\)与\(x_i=0\)相比，\(y=1\)优势比的对数值。

2. Logistic 回归建模步骤

（ 1 ）Code

利用 \(\text{scikit-learn}\) 库对这个数据建立逻辑回归模型

import pandas as pd
from sklearn.linear_model import LogisticRegression as LR
# 参数初始化
filename = '../data/bankloan.xls'
data = pd.read_excel(filename)x = data.iloc[:,:8]
y = data.iloc[:,8]lr = LR()  # 建立逻辑回归模型
lr.fit(x, y)  # 用筛选后的特征数据来训练模型
print('模型的平均准确度为：%s' % lr.score(x, y))

模型的平均准确度为：0.8085714285714286

一些参数如下：

• model.coef_ 得到[\(\beta_1\),\(\beta_2\),\(\cdots\),\(\beta_k\) ]

• model.intercept_得到\(\beta_0\),截距，默认有截距

5.1.4 决策树

决策树方法在分类、预测、规则提取等领域有着广泛应用。决策树是一种树状结构，它的每一个叶节点对应着一个分类，非叶节点对应着在某个属性上的划分，根据样本在该属性上的不同取值将其划分成若干个子集。对于非纯的叶节点，多数类的标号给出到达这个节点的样本所属的类。构造决策树的核心问题是在每一步如何选择适当的属性对样本做拆分。对一个分类问题，从已知类标记的训练样本中学习并构造出决策树是一个自上而下、分而治之的过程。

1.ID3算法

这个图片中的公式是关于信息熵的一个变体，特别是针对集合 ( A ) 中的事件 \(( a_1, a_2, \ldots, a_k )\) 对某个特定系统状态 \(S\) 的贡献。这个系统状态 \(S\) 可以进一步细分为 \(( S_1, S_2, \ldots, S_k )\)。公式中的 $E(A) $ 表示的是系统状态 $ S $ 下集合 $ A $ 中各事件的信息熵期望值。

公式 \(E(A)\) 定义如下：

\[E(A) = \sum_{j=1}^k \frac{s_{j}}{s} I(s_{j1}, s_{j2}, \ldots, s_{jm_j}) \]

这里:

• $s_j $ 表示在状态 $S_j $ 下的特定事件的数量。
• $s $是所有事件的总数。
• $I(s_{j1}, s_{j2}, \ldots, s_{jm_j}) $是在状态 $S_j $ 下，对应 $ s_j $ 事件集合的信息熵。

简而言之，这个公式计算的是给定系统状态下，各个子状态的加权信息熵总和。每个子状态的信息熵是根据该子状态的事件数量进行加权的。这种方法常用于考虑条件概率或者系统在不同子状态下的不确定性。

2. ID3 算法具体流程

示例较为简单，自行阅读。

（ 1 ）具体步骤

就是在按照具体流程一步一步计算。

（ 2 ）Code

import pandas as pd
# 参数初始化
filename = '../data/sales_data.xls'
data = pd.read_excel(filename, index_col = '序号')  # 导入数据# 数据是类别标签，要将它转换为数据
# 用1来表示“好”“是”“高”这三个属性，用-1来表示“坏”“否”“低”
data[data == '好'] = 1
data[data == '是'] = 1
data[data == '高'] = 1
data[data != 1] = -1
x = data.iloc[:,:3].astype(int)
y = data.iloc[:,3].astype(int)from sklearn.tree import DecisionTreeClassifier as DTC
dtc = DTC(criterion='entropy')  # 建立决策树模型，基于信息熵
dtc.fit(x, y)  # 训练模型# 导入相关函数，可视化决策树。
# 导出的结果是一个dot文件，需要安装Graphviz才能将它转换为pdf或png等格式。
from sklearn.tree import export_graphviz
x = pd.DataFrame(x)"""
string1 = '''
edge [fontname="NSimSun"];
node [ fontname="NSimSun" size="15,15"];
{
''' 
string2 = '}'
"""with open("../tmp/tree.dot", 'w') as f:export_graphviz(dtc, feature_names = x.columns, out_file = f)f.close()from IPython.display import Image  
from sklearn import tree
import pydotplus dot_data = tree.export_graphviz(dtc, out_file=None,  #regr_1 是对应分类器feature_names=data.columns[:3],   #对应特征的名字class_names=data.columns[3],    #对应类别的名字filled=True, rounded=True,  special_characters=True)  
graph = pydotplus.graph_from_dot_data(dot_data)  
graph.write_png('../tmp/example.png')    #保存图像
Image(graph.create_png()) 

（ 2 ）Result

5.1.5 人工神经网络

神经网络本人学习较多，所以此处不加以展开，这里用pytorch重构一下案例代码。

• 本案例数据量太少只有34条，所以使用神经网络效果肯定是不好的，这边只是实现一下代码。
• 书中给出的代码，直接在训练集上做测试，简直是不可理喻，但是此题数据量太少了，所以数据分割其实意义也不大，此处结果单纯做一个示范，不要较真。

1.Code

import pandas as pd
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix
import seaborn as sns# 参数初始化
inputfile = '../data/sales_data.xls'
data = pd.read_excel(inputfile, index_col='序号')  # 导入数据# 数据是类别标签，要将它转换为数据
# 用1来表示“好”“是”“高”这三个属性，用0来表示“坏”“否”“低”
data[data == '好'] = 1
data[data == '是'] = 1
data[data == '高'] = 1
data[data != 1] = 0# 分割数据
x_train_data = np.array(data.iloc[1:29, :3].values).astype(float)
y_train_data = np.array(data.iloc[1:29, 3].values).astype(float)
x_eval_data =  np.array(data.iloc[29:35, :3].values).astype(float)
y_eval_data =  np.array(data.iloc[29:35, 3].values).astype(float)x = torch.Tensor(x_train_data)
y = torch.Tensor(y_train_data)
x_eval = torch.Tensor(x_eval_data)
y_eval = torch.Tensor(y_eval_data)
# 创建TensorDatasets和DataLoader
dataset = TensorDataset(x, y)
dataloader = DataLoader(dataset, batch_size=10, shuffle=False)# 定义模型
class Net(nn.Module):def __init__(self):super(Net, self).__init__()self.fc1 = nn.Linear(3, 10)self.fc2 = nn.Linear(10, 1)self.relu = nn.ReLU()self.sigmoid = nn.Sigmoid()def forward(self, x):x = self.relu(self.fc1(x))x = self.sigmoid(self.fc2(x))return xmodel = Net()# 损失函数和优化器
loss_fn = nn.BCELoss()optimizer = optim.Adam(model.parameters(),lr=0.01)# 训练模型
for epoch in range(1000):for inputs, labels in dataloader:optimizer.zero_grad()outputs = model(inputs)loss = loss_fn(outputs.squeeze(), labels)loss.backward()optimizer.step()# 预测
with torch.no_grad():predictions1 = model(x_eval).squeeze()predictions = (predictions1 >= 0.5).float()  # 将输出转换为0或1# 绘制混淆矩阵
def cm_plot(original, predictions):cm = confusion_matrix(original, predictions)sns.heatmap(cm, annot=True, fmt='d', cmap='Blues')plt.xlabel('Predicted')plt.ylabel('True')plt.show()cm_plot(y_eval, predictions)

2.Result

数据量太少了，此处单纯用来练习可视化，测试数据太少不要较真。

5.1.6 分类与预测算法评价

分类与预测模型对训练集进行预测而得出的准确率并不能很好地反映预测模型未来的性能，为了有效判断一个预测模型的性能表现，需要一组没有参与预测模型建立的数据集，并在该数据集上评价预测模型的准确率，这组独立的数据集叫测试集。模型预测效果评价，通常用绝对误差与相对误差、平均绝对误差、根均方差、相对平方根误差等指标来衡量。

1. 绝对误差与相对误差

设\(Y\)表示实际值，\(\hat{Y}\)表示预测值,则称\(E\) 为绝对误差(Absolute Error)，计算公式如下：

\[E=Y-\hat{Y} \]

\(e\)为相对误差(Relative Error),计算公式如下：

\[e=\frac{Y-\hat{Y}}Y \]

有时相对误差也用百分数表示，如下式所示。

\[e=\frac{Y-\hat{Y}}Y\times100\% \]

这是一种直观的误差表示方法。

2.平均绝对误差

平均绝对误差( Mean Absolute Error, MAE),计算公式如下所示：

\[\mathrm{MAE}=\frac1n\sum_{i=1}^n\mid E_i\mid=\frac1n\sum_{i=1}^n\mid Y_i-\hat{Y_i}\mid \]

式中，\(\rm{MAE}\) 表示平均绝对误差，\(E_i\)表示第\(i\)个实际值与预测值的绝对误差，\(Y_{i}\)表示第\(i\)个实际值，\(\hat{Y}_i\)表示第\(i\)个预测值。由于预测误差有正有负，为了避免正负相抵消，故取误差的绝对值进行综合并取其平均数，这是误差分析的综合指标法之一。

3.均方误差

均方误差(Mean Squared Error, MSE),计算公式如下式所示：

\[\mathrm{MSE}=\frac1n\sum_{i=1}^nE_i^2=\frac1n\sum_{i=1}^n(Y_i-\hat{Y}_i)^2 \]

式中，\(\rm{MSE}\) 表示均方差。本方法用于还原平方失真程度。均方误差是预测误差平方之和的平均数，它避免了正负误差不能相加的问题。由于对误差\(E\)进行了平方，加强了数值大的误差在指标中的作用，从而提高了这个指标的灵敏性，是一大优点。均方误差是误差分析的综合指标法之一。

4.均方根误差

均方根误差(Root Mean Squared Error, RMSE)计算公式如下式所示。

\[\mathrm{RMSE}=\sqrt{\frac1n\sum_{i=1}^nE_i^2}=\sqrt{\frac1n\sum_{i=1}^n\left(Y_i-\hat{Y}_i\right)^2} \]

式中，\(\rm{RMSE}\) 表示均方根误差。这是均方误差的平方根，代表了预测值的离散程度，也叫标准误差，最佳拟合情况
为 \(\rm{RMSE}=0\)。均方根误差也是误差分析的综合指标之一。

5.平均绝对百分误差

平均绝对百分误差( Mean Absolute Percentage Error, MAPE)计算公式如下式所示。

\[\mathrm{MAPE}=\frac{1}{n}\sum_{i=1}^{n}\mid E_{i}\:/\:Y_{i}\mid=\frac{1}{n}\sum_{i=1}^{n}\mid(Y_{i}-\hat{Y}_{i})\:/\:Y_{i}\mid \]

式中，\(\rm{MAPE}\) 表示平均绝对百分误差。一般认为 MAPE 小于 10 时，预测精度较高。

6. Kappa 统计

7. 识别准确度

8. 识别精确率

9. 反馈率

10. ROC 曲线

11. 混淆矩阵

神经网络那个例子里的结果就是混淆矩阵，看一下就明白了。

5.1.7 Python 分类预测模型特点