经过数据探索与数据预处理,得到了可以直接建模的数据,根据挖掘目标和数据形式可以建立分类与预测,聚类分析,关联规则,时序模式和偏差检测等模型,帮助企业提取数据中的商业价值.
分类与预测
对于之前的背景假设,餐饮企业会经常遇到如下问题:
- 如何基于菜品历史销售情况,以及节假日,气候,竞争对手等影响因素,对菜品销量进行趋势预测?
- 如何预测未来一段时间哪些顾客会流失,哪些顾客会最有可能成为vip?
- 如何预测一种新产品的销售量,以及在那种类型中的顾客会较受欢迎?
除此之外,餐厅经理需要通过数据分析知道具有某些特征的顾客的消费习惯,老板希望知道下一个月的销售收入,原材料需要投入,这些都是分类与预测的例子.
实现过程
- 分类是构造一个分类模型,输入样本的属性值,输出对应类别,将每一个样本映射到预先定义好的类别.分类模型建立在已有类标记的数据集上,模型在已有样本的准确率可以方便计算,属于有监督的学习.
- 预测是指建立两种或两种以上变量间的相互依赖的函数模型,然后进行预测或控制,
以上两种实现过程类似,以分类模型为例,分类算法有两步骤:第一步是学习步,通过归纳分析样本集建立分类模型得到分类规则;第二步是分类步,先用已知的测试样本集评估分类规则的准确率,如果准确率是可以接受的,则使用该模型对未知类标号的待测样本集进行预测.预测模型实现也有两步:第一步通过训练集建立预测属性(数值型)的函数模型,第二步在模型通过检验后进行预测或控制,如图所示:
常用分类与预测算法
常用的分类与预测算法简介如下表:
| 算法名称 | 算法描述 |
|---|---|
| 回归分析 | 回归分析是确定预测属性(数值型)与其他变量间互相依赖的定量关系最常用的统计学方法,包括线性回归,非线性回归,Logistic回归,岭回归,主成分回归,偏最小二乘回归等. |
| 决策树 | 决策树采用自顶向下的递归方式,在内部节点进行预测属性值的比较,并根据不同的属性值从该节点向下分支最终得到划分类的叶节点 |
| 人工神经网络 | 人工神经网络是一种模仿大脑神经网络结构和功能而建立的信息处理系统,表示神经网络的输入与输出变量之间关系的模型 |
| 贝叶斯网络 | 又称信度网络,是Bayes方法的扩展,是目前不确定知识表达和推理领域最有效的理论模型 |
| 支持向量机 | 通过某种非线性映射,把低维非线性转化为高维线性可分,在高维空间进行线性分析的算法 |
回归分析
回归分析是通过建立模型来研究变量之间相互关系的密切程度,结构状态以及进行模型预测的一种有效工具,回归分析研究范围大致如下:
主要回归模型分类:
| 回归模型名称 | 适用条件 | 算法描述 |
|---|---|---|
| 线性回归 | 因变量与自变量是线性关系 | 对一个或多个自变量与因变量之间的线性关系进行建模,可用最小二乘法求解模型系数 |
| 非线性回归 | 因变量与自变量不都是线性关系 | 对一个或多个自变量与因变量之间的非线性关系进行建模,如果非线性相关可以通过简单的函数变化转化成线性相关,则用线性相关求解,否则使用非线性最小二乘法求解 |
| Logistic回归 | 因变量一般有1和0两种取值 | 广义线性回归的特例,利用Logistic函数将因变量的取值范围控制在0和1之间,表示取值为1的概率 |
| 岭回归 | 参与建模的自变量之间具有多重共线性 | 改进的最小二乘估计的方法 |
| 主成分回归 | 参与建模的自变量之间具有多重共线性 | 根据主成分分析思想提出来的,是对最小二乘法的改进,参数估计的一种有偏估计,可消除自变量间的多重共线性 |
线性回归是相对简单的回归模型,但是因变量与自变量之间呈现某种曲线关系,就需要建立非线性回归模型.Logistic回归属于概率型非线性回归,分为二分类和多分类的回归模型,对于二分类的Logistic回归,因变量y只有"是,否"两个取值,记为1和0.假设在自变量X1,X2,…,Xp,作用下,y取"是"的概率是p,取否概率是(1-p),研究的是当y取"是"发生概率p与自变量X1,X2,…,Xp的关系.当自变量之间出现多重共线性时,用最小二乘估计回归系数将不会准确,消除多重共线性参数改进的估计方法主要有岭回归和主成分回归.下面常用二分类Logistic回归模型原理进行解释.
- Logistic回归分析介绍
Logistic函数:Logistic回归模型中因变量的只有1-0(如是和否,发生和不发生)两种取值,假设在p个独立自变量X1,X2,…,作用下,记y取1的概率是p=P(y=1|x),取0的概率是1-p,取1和0的概率之比为p/1-p,称为事件的优势比(odds),对优势比取自然对数即取得Logistic变换ln(p/1-p),令Logit§=ln(p/1-p)=z,则Logistic函数为:
当p在(0,1)之间变化时,odds的取值范围是0到正无穷,则ln(p/1-p)de取值范围是实数.
Logistic回归模型:是建立ln(p/1-p)与自变量的线性回归模型. Logistic回归模型为:
Logistic回归模型解释,β0即为没有自变量,即X1,…,Xp全部取0,为概率之比的自然对数,βi为自变量变化时,Xi=1与Xi=0相比,y=1的优势比对数值. - Logistic回归建模步骤
根据分析目标设置指标变量(因变量和自变量),然后收集数据,根据收集到的数据,对特征再次筛选;y取1的概率是p,y取0的概率是1-p,列出线性回归方程,估计模型的回归系数;进行模型检验,根据正确率,其次有混淆矩阵,ROC曲线,KS值等;模型应用,输入自变量获得预测值,或通过预测变量控制自变量值,具体如图:
下面是某银行在降低贷款拖欠率的数据进行逻辑回归建模,该数据示例在demo/data/bankloan.xls,利用Scikit-Learn对这个数据进行逻辑回归分析,首先进行特征筛选,特征筛选的方法有很多,主要包含在feature_selection库中,比较简单的是通过F检验(f_regression)来给出各个特征的F值和P值,从而可以筛选变量(选择F值大的或者p值小的).其次有递归特征消除(RFE)和稳定性选择等比较新的方法,然后利用筛选的特征(提高准确性)建立逻辑回归模型,其代码:
#-*- coding: utf-8 -*-
#逻辑回归 自动建模
import pandas as pd
#参数初始化
filename = '../data/bankloan.xls'
data = pd.read_excel(filename)
x = data.iloc[:,:8].as_matrix()#读取自变量,并处理为pd的二维数组，前是行，下标全部的行；，后是列，列选下标0-7
y = data.iloc[:,8].as_matrix()#读取因变量,是否违约,行全选，列选下标8
#使用稳定性选择方法中的随机逻辑回归进行特征筛选，利用筛选后的特征建立逻辑回归模型，输出平均正确率
from sklearn.linear_model import LogisticRegression as LR
lr = LR() #建立Logistic回归分析模型
lr.fit(x, y) #用筛选后的特征数据来训练模型
print(u'逻辑回归模型训练结束。')
print(u'模型的平均正确率为：%s' % lr.score(x, y)) #给出模型的平均正确率，本例为80.57%
因为Scikit-Learn库版本更新,并没有随机逻辑回归函数,所以上例并没有采用随机逻辑回归剔除变量,所以模型正确率低了1%左右,使用逻辑回归筛选出来的变量与结果具有比较强的线性相关性,然而被筛选掉的变量也可能具有非线性关系,因此,还要根据背景对筛选结果进行分析,对于非线性关系,筛选方法有决策树,神经网络等.
决策树
决策树方法在分类,预测,规则提取等领域有着广泛应用,决策树是一树状结构,每一个叶节点对应一个分类,非叶节点对应着某个属性的划分,根据样本在该属性的不同取值划分成若干子集.构造决策树的核心问题是每一步如何选择适当的属性对样本做拆分,对一个分类问题,从已知标记的分类样本中构造出决策树是一个自上而下,分而治之的过程,常用决策树方法有:
| 决策树算法 | 算法描述 |
|---|---|
| ID3算法 | 其核心是在决策树的各节点上,使用信息增益的方法作为属性的选择标准,来帮助确定生成的每一个节点所应采用的合适属性 |
| C4.5算法 | C4.5决策树生成算法相对于ID3算法改进是使用信息增益率来选择节点属性,ID3只可以处理离散的描述属性,而C4.5不仅可以处理离散的描述属性,还可以处理连续的描述属性 |
| CART算法 | CART决策树是有效的非参数分类和回归方法,通过构建树,修剪树,评估树来构建二叉树,当终结点是连续变量时,称为回归树,当终结点是分类变量,该树为分叉树 |
-
ID3算法简介及基本原理
熵定义为信息的期望值,在信息论与概率论统计中，熵表示随机变量不确定性的度量.
ID3算法基于信息熵来选择最佳测试属性,它选择当前样本集中具有最大信息增益值的属性作为测试属性,样本集的划分基于测试属性的取值进行,测试属性拆分有多少不同取值就将样本集划分为多少子样本集,同时决策树上相应于该样本集的节点长出新的叶子节点.ID3算法用信息增益值度量不确定性,增益越大,不确定越小,因此.每一个非叶节点选择信息增益最大的属性可以得到最确定的分类,得到较小的决策树.
设S是s个抽样数据样本的集合,假定类别属性有m个不同的值;Ci(i=1,2,…,m),也就是说S中C属性有m个分类.设si是类别的Ci的样本数,某一分类个数,其中Pi是任意样本属于Ci的概率,一般通过si/s估计它的该类总信息熵为:
设一个新的属性A具有k个不同的值Ai(i=1,2,…,k),利用属性A将集合S划分为k个子集(S1,S2,…,Sk),其中Sj包含集合S中属性A取Aj的样本.若选择A作为测试属性,则子集作为从集合S的节点生长出来的叶子节点.其中将Sj继续进行分类,则基于C属性的话,Sij是子集Sj中Ci的样本集,则根据属性A划分类别样本总信息熵值为:
最后利用属性A划分样本集S的信息增益为:
说明选择属性A的信息熵越大,信息增益越小,不确定性越小,对应的属性A的子集通过递归调用,生成其他属性作为子节点和分支来生成整个决策树. -
ID3算法具体流程
ID3算法具体实现步骤如下:对当前样本的集合,计算所有属性的信息增益;选择信息增益最大的属性作为测试属性,把测试属性取值相同的样本划分为同一个子集样本;若单个子集样本的类别属性只有单个属性,则分支为叶子节点,判断其属性值并标上相应符号,否则对样本子集继续调用本算法.
结合餐饮案例了解ID3的具体步骤,我们预测销量需要构建模型分析天气,是否周末和是否有促销活动对销量的影响,下面以单个门店简单为例,将不下雨的天气设置为好,下雨雪的为坏,周末设置为是,非周末设置为否,促销有无设置为是否,将销量数据离散化,大于均值为高,低于均值为低,经过处理得到数据demo/data/sales_data.xls.
采用ID3算法构建决策树模型步骤为:
获得分类样本高低总信息熵值,数据总记录34,高为18,低有16:
计算每一个测试属性的信息熵,以天气为例,其属性值有好坏之分,其中天气为好条件下,高的销量为11,低为6,表示(11,6);天气为坏,销量为高为7,低为10,可表示为(7,10),则天气信息熵计算过程:
同理计算周末,促销和销量的信息熵,再计算天气,周末,促销属性的信息增益值:
可以得出周末属性的信息增益最大,它的两个属性值是否作为该根节点的两个分支,然后对每一个分支进行信息增益计算,如此循环往复,直到没有新节点分支,最终构成的决策树为:
可以从决策树提炼规则,例如周末好天气销售高,不是周末坏天气促销销量低之类的.
ID3决策树算法会偏向选择取值较多的,所谓高度分支属性,这一类属性不一定是最优属性,此外只能处理离散属性,针对连续属性需要离散化,为了解决这些问题,人们采用信息增益率C4.5算法,此外常用决策树算法还有CART算法,SLIQ算法,SPRINT算法和PUBLIC算法等,使用Sickit-Learn建立基于信息熵的决策树模型代码
#-*- coding: utf-8 -*-
#使用ID3决策树算法预测销量高低
import pandas as pd
#参数初始化
inputfile = '../data/sales_data.xls'
data = pd.read_excel(inputfile, index_col = u'序号') #导入数据
#数据是类别标签，要将它转换为数据
#用1来表示“好”、“是”、“高”这三个属性，用-1来表示“坏”、“否”、“低”
data[data == u'好'] = 1
data[data == u'是'] = 1
data[data == u'高'] = 1
data[data != 1] = -1
x = data.iloc[:,:3].as_matrix().astype(int)
y = data.iloc[:,3].as_matrix().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)
from sklearn.externals.six import StringIO
x = pd.DataFrame(x)
with open("tree.dot", 'w') as f:
f = export_graphviz(dtc, feature_names = x.columns, out_file = f)
运行代码后生成tree.dot文件,添加一下两行
digraph Tree {
edge [fontname="SimHei"];
node [fontname="SimHei"] ;
0 [label="1 <= 0.0\nentropy = 0.998\nsamples = 34\nvalue = [16, 18]"] ;
保存为UTF-8格式,安装Graphviz,然后在命令行编译命令:
dot -Tpdf tree.out -o tree.pdf
dot -Tpng tree.out -o tree.pdf
生成结果如图:
人工神经网络
人工神经网络是模拟生物神经网络进行信息处理的一种数学模型,人工神经元是人工神经网络操作处理的基本信息处理单位,它是人工神经网络设计基础.一个人工神经元对输入信号X={X1,X2,…,Xn}依据其权重W={W1,W2,…,Wn}求出输入为{W1X1+W2X2,…,WnXn},然后然后被一个数学函数用来计算决定是否激发神经元,还有一个函数（也许是不变，就是复制）计算人工神经元的输出（有时依赖于某个门限）,人工神经网络把这些人工神经元融合一起用于处理信息,如图所示:
人工神经网络的学习也称为训练,指的是神经网络在外部环境刺激下调整神经网络参数,使神经网络以一种新的方式对外部环境做出反应的过程.在分类与预测中,根据指定的训练样本,调整人工神经网络参数以使网络输出接近已知的样本类标记或其他形式的因变量,激活函数分类表如下:
在人工神经网络的发展过程中,提出了多种不同的学习算法,没有一种特定的学习算法适用于所有的网络结构,在分类与预测中,δ学习规则(误差校正学习算法)是使用最广泛的一种.误差校正学习算法根据神经元的输出误差对神经元的连接强度进行修正,属于有指导的学习,设神经元i为输入,神经元j为输出,他们的连接权值为Wij,则对权值修正为ΔWij=ηδjYi,其中η为学习率,δj=Tj-Yj,为j的偏差,即输出神经元j的期望输出和实际输出之差,示意图:
神经网络训练是否完成常用误差函数(目标函数)E来衡量,当误差函数小于某个设定的值即停止神经网络的训练,误差函数为衡量实际输出向量Yk与期望值向量Tk误差大小的函数,常用二乘误差函数来定义:
使用人工神经网络模型需要确定网络连接的拓扑结构,神经元的特征,学习规则等,目前,已有40种人工神经网络模型,常用实现分类和预测的人工神经网络算法见表:
| 算法名称 | 算法描述 |
|---|---|
| BP神经网络 | 是一种按照误差逆传播算法训练的多层前馈网络,学习算法是δ学习规则,是目前应用最广泛的神经网络模型之一 |
| LM神经网络 | 是基于梯度下降法和牛顿法结合的多层前馈网络,特点是:迭代次数少,收敛速度快,精确度高 |
| RBF径向基神经网络 | RBF神经网络能够以任意精度逼近任意连续函数任意连续函数,从输入层到隐含层的变换是非线性的,而从隐含层到输出层的变换是线性的,适合解决分类问题 |
| FNN模糊神经网络 | FNN模糊神经网络是具有模糊系数或者输入信号是模糊量神经网络,是模糊系统和神经网络相结合的产物,它汇聚了神经网络和模糊系统的优点,集联想,识别,自适应以及模糊信息处理集一体 |
| GMDH神经网络 | GMDH神经网络也称为多项式网络,他是前馈神经网络用于预测的网络,它的特点是网络结构不固定,而且在训练过程中不断改变 |
| ANFIS自适应神经网络 | 在一个全部模糊的结构中,不知不觉中向训练集学习,自动产生,修正,并高度概括最佳输入与输出变量的隶属函数及模糊规则 |
反向传播算法(BP)的特征是利用输出后的误差来估计输出层的直接前导层的误差,再用这个误差估计更前一层的误差,如此一层一层反向传播下去就获得了所有其他各层的误差估计,这就形成了以输出层表现出的误差沿着与输入传送相反的方向逐级向网络输入层传递的过程,以典型的三层BP网络为例,描述标准BP算法,图示有三个输入节点,4个隐层节点,1个输出节点的三层BP神经网络.
BP算法的学习过程由信号的正向传播与误差的逆向传播两个过程组成,正向传播时,输入信号经隐层处理后,传向输出层.若输出层节点未能得到期望的输出,则转入误差逆向传播阶段,将输出误差按照某种子形式,通过隐层向输入层返回,并分摊给隐层的4个节点与输入层的3个输入节点,从而获得各层单元的参考误差(误差信号),作为修改各单元权值的依据,这种信号的正向传播与误差的逆向传播的各层权矩阵的修改过程,也就是网络学习(训练)的过程,在此过程一直进行到网络输出的误差逐渐减少到可接受的程度或者设定的学习次数为止学习过程流程图为:
算法开始后,给定学习次数上线,初始化学习次数为0,对权值和阈值赋值较小的随机数,一般在[-1,1]之间,输入样本数据,网络正向传播,得到中间层与输出层的值,比较输出层与教师信号值的误差,通过误差函数判断是否小于误差上限,如不小于误差上限,对中间层和输出层的权值和阈值进行更新,更新算法为δ学习规则,在次将样本数据作为输入,得到中间层和输出层的值,计算误差E是否小于上限,学习次数是否达到指标值,如果达到,学习结束.
BP算法只用到了均方误差函数对权值阈值的一阶导数信息,使算法存在收敛速度慢,易陷入局部极小等缺陷,为了解决这一问题,2006年提出非监督贪心逐层训练算法,并以此基础发展为深度学习算法.Scikit-Learn并没有提供神经网络模型,所以我们使用Keras,这是深度学习库,我们建立神经网络模型有3个输入节点,10个隐藏节点和1个输出节点,其代码:
#-*- coding: utf-8 -*-
#使用神经网络算法预测销量高低
import pandas as pd
#参数初始化
inputfile = '../data/sales_data.xls'
data = pd.read_excel(inputfile, index_col = u'序号') #导入数据
#数据是类别标签，要将它转换为数据
#用1来表示“好”、“是”、“高”这三个属性，用0来表示“坏”、“否”、“低”
data[data == u'好'] = 1
data[data == u'是'] = 1
data[data == u'高'] = 1
data[data != 1] = 0
x = data.iloc[:,:3].as_matrix().astype(int)
y = data.iloc[:,3].as_matrix().astype(int)
from keras.models import Sequential
from keras.layers.core import Dense, Activation
model = Sequential() #建立模型
model.add(Dense(input_dim = 3, output_dim = 10))
model.add(Activation('relu')) #用relu函数作为激活函数，能够大幅提供准确度
model.add(Dense(input_dim = 10, output_dim = 1))
model.add(Activation('sigmoid')) #由于是0-1输出，用sigmoid函数作为激活函数
model.compile(loss = 'binary_crossentropy', optimizer = 'adam')
#编译模型。由于我们做的是二元分类，所以我们指定损失函数为binary_crossentropy，以及模式为binary
#另外常见的损失函数还有mean_squared_error、categorical_crossentropy等，请阅读帮助文件。
#求解方法我们指定用adam，还有sgd、rmsprop等可选
model.fit(x, y, nb_epoch = 1000, batch_size = 10) #训练模型，学习一千次
yp = model.predict_classes(x).reshape(len(y)) #分类预测
from cm_plot import * #导入自行编写的混淆矩阵可视化函数
cm_plot(y,yp).show() #显示混淆矩阵可视化结果
运行以上代码,可以看出检测34个样本,预测正确的个数为26个,预测正确率为76.4%,正确率低的原因是训练需要较多的样本,由于训练样本缺少造成,神经网络容易出现过拟合现象,目前深度学习网络通过随机的让部分神经网络节点休眠防止过拟合.
分类与预测算法评价
分类与预测模型对训练集进行预测而得出的准确率并不能很好地反映预测模型的性能,所以为了有效判断预测模型性能表现,我们需要一组没有参与预测模型建立的数据集,并在该数据集上评价预测模型的准确率,这组独立的数据集叫测试集.模型预测效果评价,通常用相对/绝对误差,平均绝对误差,均方误差,均方根误差等指标衡量.
- 绝对误差与相对误差,实际值与预测值的差叫做绝对误差(E),计算公式:
相对误差e计算公式:
也可以用百分比表示:
- 平均绝对误差,计算公式如下,为避免正负相抵消,故取误差的绝对值进行综合并取其平均数,这是误差分析的综合指标法之一:
- 均方误差,用于还原平方失真程度,对误差E进行平方,加强了数值大的误差在指标中的作用,从而提高了指标的灵敏性,是一大优点,均方误差是误差分析的综合指标法之一.
- 均方根误差,代表了预测值的离散程度,也称标准误差,最佳拟合情况RMSE=0,均方根误差也是误差分析的综合指标之一.
- 平均绝对误差,MAPE表示平均绝对百分误差,一般认为MAPE小于10时,预测精度较高.
- 混淆矩阵,混淆矩阵（Confusion Matrix），它的本质远没有它的名字听上去那么拉风。矩阵，可以理解为就是一张表格，混淆矩阵其实就是一张表格而已。 以分类模型中最简单的二分类为例，对于这种问题，我们的模型最终需要判断样本的结果是0还是1，或者说是positive还是negative。我们通过样本的采集，能够直接知道真实情况下，哪些数据结果是positive，哪些结果是negative。同时，我们通过用样本数据跑出分类型模型的结果，也可以知道模型认为这些数据哪些是positive，哪些是negative。因此，我们就能得到这样四个基础指标，我称他们是一级指标（最底层的）：
真实值是positive，模型认为是positive的数量（True Positive=TP）
真实值是positive，模型认为是negative的数量（False Negative=FN）：这就是统计学上的第一类错误（Type I Error）
真实值是negative，模型认为是positive的数量（False Positive=FP）：这就是统计学上的第二类错误（Type II Error）
真实值是negative，模型认为是negative的数量（True Negative=TN）
将这四个指标一起呈现在表格中，就能得到如下这样一个矩阵，我们称它为混淆矩阵（Confusion Matrix）：
- Kappa统计,Kappa统计是比较两个或多个观测者对同一事物,或观测者对同一事物的两次或多次的观测结果是否一致,以由于机遇造成的一致性和实际观测值一致性之间的差别大小作为评价基础的统计指标.Kappa统计量和加权Kappa统计量不仅可以用于无序和有序分类变量资料的一致性,重现性检验,而且可以给出一个反映一致性大小的"量"值.Kappa取值在[-1,+1]之间,其值的大小均有不同意义.
它是通过把所有地表真实分类中的像元总数（N）乘以混淆矩阵对角线（Xkk）的和，再减去某一类地表真实像元总数与被误分成该类像元总数之积对所有类别求和的结果，再除以总像元数的平方减去某一类中地表真实像元总数与该类中被误分成该类像元总数之积对所有类别求和的结果所得到的。
——来自百科
kappa计算结果为[-1,1],但通常kappa是落在 0-1 间,可分为五组来表示不同级别的一致性:0.0-0.20极低的一致性(slight)、0.21-0.40一般的一致性(fair)、0.41-0.60 中等的一致性(moderate)、0.61-0.80 高度的一致性(substantial)和0.81-1几乎完全一致(almost perfect).计算公式：
po是每一类正确分类的样本数量之和除以总样本数，也就是总体分类精度,假设每一类的真实样本个数分别为a1,a2,…,aC,而预测出来的每一类的样本个数分别为b1,b2,…,bC,总样本个数为n,则有：pe=a1×b1+a2×b2+…+aC×bC / n×n.运算举例,为了更好的理解上述运算的过程，这里举例说明一下：学生考试的作文成绩，由两个老师给出 好、中、差三档的打分，现在已知两位老师的打分结果，需要计算两位老师打分之间的相关性kappa系数：
上图其实是混淆矩阵,可以理解对角线为评分一致,其他为A老师中,B老师为好之类不一致
Po = (10+35+15) / 87 = 0.689
a1 = 10+2+8 = 20; a2 = 5+35+5 = 45; a3 = 5+2+15 = 22;
b1 = 10+5+5 = 20; b2 = 2+35+2 = 39; b3 = 8+5+15 = 28;
Pe = (a1*b1 + a2*b2 + a3*b3) / (87*87) = 0.455
K = (Po-Pe) / (1-Pe) = 0.4293578
这样我们就得到了kappa系数。
Kappa = +1说明两次判断结果完全一致.
Kappa = -1说明两次判断结果完全不一致.
Kappa = 0说明两次判断的结果是机遇造成.
Kappa < 0说明一致程度比机遇造成的还差,两次检测结果不一致,在实际应用中无意义.
Kappa > 0说明有意义,Kappa越大,说明一致性越好.
Kappa >= 0.75说明已经取得相当满意的一致程度.
Kappa < 0.4 说明一致程度不够
8. 识别准确度,识别准确度定义如下,
8. 识别精确率,定义如下,
10. 反馈率,定义如下,
11. ROC曲线,受试者工作特征曲线是一种非常有效的模型评价方法,在机器学习领域,ROC曲线的横坐标为假阳性率(FPR),FPR = FP/N,纵坐标为真阳性率(TPR),TPR=TP/P,P是真实的正样本的数量,N是真实的负样本的数量,TP是P个正样本中被分类器预测为正样本的个数,FP是N个负样本中被分类器预测为正样本的个数,假如雷达显示器上出现了10次信号,其中3次确实是敌机来袭(P=3),另外7次都是飞鸟的信号(N=7),雷达兵判断3次敌机来袭中有2次是敌机来袭(TP=2),那么真阳性率TPR=TP/P=2/3,对于7次飞鸟信号,有1次误判成敌机来袭(FP=1),那么假阳性率FPR=FP/N=1/7,所以对于这个雷达兵分类器来说,这组分类结果应该是ROC曲线上的点(1/7,2/3).ROC曲线下的面积叫AUC,AUC的值能够量化的反映基于ROC曲线衡量出的模型性能.AUC等于随机挑选一个正样本和负样本是,分类器将正样本排在前面的概率,AUG其值接近1,说明算法效果越好.
Python分类预测模型特点
Python程序是面向对象的,放到建模之中,不管是在Scikit-Learn还是在Keras中,建模第一步是建立一个对象,这个对象是空白的,需要进一步训练,然后我们设置模型参数,通过fit方法对模型进行训练,最后通过predict方法预测结果,当然还有一些对完成模型的评估,如score等,常见的模型评价和在Python中实现:
| 模型 | 模型特点 | 位于 |
|---|---|---|
| 逻辑回归 | 比较基础的线性分类模型 | sklearn.linear_model |
| SVM | 强大的模型,可以用来回归,预测,分类等,而根据选择不同的核函数,模型可以是线性/非线性 | sklearn.svm |
| 决策树 | 基于分类讨论,逐步细化的分类模型 | sklearn.tree |
| 随机森林 | 与决策树类似,精度高于决策树,缺点是随机性丧失决策树可解释性 | sklearn.ensemble |
| 朴素贝叶斯 | 基于概率的思想简单有效的分类模型,给出容易理解的概率解释 | sklearn.naive_bayes |
| 神经网络 | 具有强大的拟合能力,用于拟合分类,是深度学习基础 | Keras |
聚类分析
餐饮企业会面临这样的问题,如何通过餐饮客户消费行为的测量,进一步评判客户的价值和对客户进行细分,找到有价值的客户和需要关注的用户,对菜品进行分析,以便区分哪些菜品畅销毛利高,哪些菜品毛利低,餐饮企业遇到这些问题,可以通过聚类分析解决.
常用聚类分析算法
聚类分析是在没有给定划分类别的情况下,根据数据相似度进行样本分组的一种方法,与分类模型需要使用有类标记样本构成的训练数据不同,聚类模型可以建立在无类标记的数据上,是一种非监督的算法模型.聚类分类根据数据自身的距离或者相似度划分为若干组,划分的原则是组内距离最小化而组间(外部)距离最大化,如图所示,
常用聚类方法:
| 类别 | 包括的主要算法 |
|---|---|
| 划分(分裂)方法 | K-Means算法(K-平均),K-MEDOIDS算法(K-中心点),CLARANS算法(基于选择的算法) |
| 层次分析方法 | BIRCH算法(平衡迭代规约和聚类),CURE算法(代表点聚类),CHAMELEON算法(动态模型) |
| 基于密度方法 | DBSCAN算法(基于高密度连续区域),DENCLUE算法(密度分布函数),OPTICS算法(对象排序识别) |
| 基于网络算法 | STING算法(统计信息网络),CLIOUE(聚类高维空间),WAVE-CLUSTER算法(小波变换) |
| 基于模型的方法 | 统计学方法,神经网络方法 |
常用聚类算法表:
| 算法名称 | 算法描述 |
|---|---|
| K-Means | K-均值聚类也称快速聚类法,在最小化误差函数的基础上将数据划分为预定的类数K,该算法原理简单并且易于处理大量数据 |
| K-中心点 | K-均值算法对孤立的点的敏感性,K-中心点不采用簇中对象的平均值作为簇中心,而选用簇中离平均值最近的对象作为簇中心 |
| 系统聚类 | 系统聚类也称多层次聚类,分类的单位由高到低呈树形结构,且所处结构位置越低,其所包含的对象就越少,但这些对象的共同特征越多,该聚类方法只适合在小数据量的时候使用,数据量大的话非常慢 |
K-Means聚类算法
K-Means算法是典型的基于距离的非层次聚类算法,在最小化误差函数的基础上将数划分为预定的类数K,采用距离作为相似性的评价指标,即认为两个对象距离越近,相似度越大.
- 算法过程
1.从N个样本数据中随机选择K个对象作为初始的聚类中心.
2.分别计算每一个样本到各个聚类中心的距离,将对象分配到距离最近的聚类中.
3.所有对象分配完毕后,重新计算K个聚类的中心.
4.与前一次计算得到的K个聚类中心比较,如果聚类中心发生变化,转过程2,否则过程5
5.当质心不发生变化时停止并输出聚类结果.
聚类的选择可能依赖于初始聚类中心的随机选择,可能使得结果偏离全局最优分类,在实践中,为了得到较好的结果,通常选择不同的聚类中心,多次运行算法.在所有对象分配完成,重新计算K个聚类中心时,对于连续数据聚类中心取该簇的的均值,当样本的某些属性是分类变量时,均值可能无定义,可使用K-众数方法. - 数据类型与相似性度量
1.连续属性,对于连续属性,要对各属性值进行零-均值规范,在进行距离的计算,在K-Means聚类算法中,一般需要度量样本之间的距离,样本与簇之间的距离以及簇与簇之间的距离.
度量样本间的相似性最常用的是欧几里得距离,曼哈顿距离和闵科夫斯基距离,样本与簇之间的距离可以用样本样本到簇中心的距离d(ei,x),簇与簇之间的距离可以用簇中心的距离d(ei,ej).
用p个属性来表示n个样本的数据矩阵为:
欧几里得距离:
曼哈顿距离:
闵科夫斯基距离,q为正整数,q=1为曼哈顿距离;q=2时为欧几里得距离:
2.文档数据,对于文档数据使用余弦相似性度量,现将文档数据整理成文档-词矩阵格式:
两个文档间的相似度计算公式为:
- 目标函数
使用误差平方和和SEE作为度量聚类质量的目标函数,对于两种不同的聚类结果,选择误差平方和较小的分类结果.
连续属性的SSE计算公式:
文档数据SSE计算公式:
簇Ei的聚类中心ei计算公式为:
上述公式的符号说明:
| 符号 | 含义 |
|---|---|
| K | 聚类簇个数 |
| Ei | 第i个簇 |
| x | 对象(样本) |
| ei | 簇Ei的聚类中心 |
| ni | 第i个簇中样本个数 |
下面结合餐饮企业分析,具体消费行为特征数据在demo/data/consumption_data.xls,ID为用户ID,R为最近一次消费时间间隔,F为消费频率,M为消费总金额,根据这些数据将客户分类成不同客户群,并评价客户群价值.我们采用K-Means聚类算法,设定聚类个数K为3,最大迭代次数为500次,距离函数采取欧式距离,代码:
import pandas as pd
#参数初始化
inputfile = '../data/consumption_data.xls' #销量及其他属性数据
outputfile = '../tmp/data_type.xls' #保存结果的文件名
k = 3 #聚类的类别
iteration = 500 #聚类最大循环次数
data = pd.read_excel(inputfile, index_col = 'Id') #读取数据
data_zs = 1.0*(data - data.mean())/data.std() #数据标准化
from sklearn.cluster import KMeans
model = KMeans(n_clusters = k, n_jobs = 4, max_iter = iteration) #分为k类，并发数4
model.fit(data_zs) #开始聚类
#简单打印结果
r1 = pd.Series(model.labels_).value_counts() #统计各个类别的数目
r2 = pd.DataFrame(model.cluster_centers_) #找出聚类中心
r = pd.concat([r2, r1], axis = 1) #横向连接（0是纵向），得到聚类中心对应的类别下的数目
r.columns = list(data.columns) + [u'类别数目'] #重命名表头
print(r)
#详细输出原始数据及其类别
r = pd.concat([data, pd.Series(model.labels_, index = data.index)], axis = 1) #详细输出每个样本对应的类别
r.columns = list(data.columns) + [u'聚类类别'] #重命名表头
r.to_excel(outputfile) #保存结果
def density_plot(data): #自定义作图函数
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号
p = data.plot(kind='kde', linewidth = 2, subplots = True, sharex = False)
[p[i].set_ylabel(u'密度') for i in range(k)]
plt.legend()
return plt
pic_output = '../tmp/pd_' #概率密度图文件名前缀
for i in range(k):
density_plot(data[r[u'聚类类别']==i]).savefig(u'%s%s.png' %(pic_output, i))
输出结果:
聚类算法输出结果三类:
R F M 类别数目
0 -0.149353 -0.658893 -0.271780 559
1 -0.160451 1.114802 0.392844 341
2 3.455055 -0.295654 0.449123 40
获得分群一概率密度图:
分群一客户价值分析,消费间隔相对较小,主要集中在0-30天,消费次数集F中在10-25次,消费金额在500到2000.
获得分群二概率密度图:
分群二客户价值分析,消费间隔分布在0-30天,消费次数集F中在0-12次,消费金额在0到1800.
获得分群三概率密度图:
分群三客户价值分析,消费间隔分布在30-80天,消费次数集F中在0-15次,消费金额在0到2000.
对比分析,分群一时间消费时间间隔短,次数多,而且消费金额大,是高消费,高价值人群,分群二是一般客户,分群三是价值较低的客户群体.
聚类算法评价
聚类分析近根据样本数据本身将样本数据进行分组,其目标是为了实现组内对象相互之间是相似的(相关的),而不同组中的对象是不同的(不相关).组内相似性越大,组间差别越大,聚类效果越好.
具体聚类效果评估方法列举介绍(具体查看周志华机器学习)：
- purity,purity方法是极为简单的一种聚类评价方法，只需计算正确聚类的文档数占总文档数的比例： 其中Ω = {ω1,ω2, . . . ,ωK}是聚类的集合ωK表示第k个聚类的集合。C = {c1, c2, . . . , cJ}是文档集合，cJ表示第J个文档。N表示文档总数。
如上图的purity = ( 3+ 4 + 5) / 17 = 0.71,其中第一类正确的有5个，第二个4个，第三个3个，总文档数17。purity方法的优势是方便计算，值在0～1之间，完全错误的聚类方法值为0，完全正确的方法值为1。同时，purity方法的缺点也很明显它无法对退化的聚类方法给出正确的评价，设想如果聚类算法把每篇文档单独聚成一类，那么算法认为所有文档都被正确分类，那么purity值为1！而这显然不是想要的结果。 - RI,实际上这是一种用排列组合原理来对聚类进行评价的手段，公式如下：
其中TP是指被聚在一类的两个文档被正确分类了，TN是只不应该被聚在一类的两个文档被正确分开了，FP只不应该放在一类的文档被错误的放在了一类，FN只不应该分开的文档被错误的分开了。对上图
TP＋FP ＝ C(2,6) + C(2,6) + C(2,5) = 15 + 15 + 10 = 40 其中C(n,m)是指在m中任选n个的组合数。
TP = C(2,5) + C(2,4) + C(2,3) + C(2,2) = 20
FP = 40 - 20 = 20
相似的方法可以计算出TN = 72 FN = 24
所以RI ＝ ( 20 + 72) / ( 20 + 20 + 72 +24) = 0.68 - F值,这是基于上述RI方法衍生出的一个方法，
RI方法有个特点就是把准确率和召回率看得同等重要，事实上有时候我们可能需要某一特性更多一点，这时候就适合F值方法.
Python主要聚类分析算法
Python聚类相关算法集中在Scikit-Learn中,Python里面实现聚类主要包括K-Means聚类,层次聚类,FCM以及神经网络聚类,其主要函数列表:
| 对象名 | 函数功能 | 所属工具箱 |
|---|---|---|
| KMeans | K均值聚类 | sklearn.cluster |
| AffinityPropagation | 吸引力传播聚类,几乎优于所有其他方法,不需要指定聚类数,但运行效率较低 | sklearn.cluster |
| MeanShift | 均值漂移聚类算法 | sklearn.cluster |
| SpectralClustering | 谱聚类,具有效果比K均值好,速度比K均值快 | sklearn.cluster |
| AgglomerativeClustering | 层次聚类,给出一棵聚类层次树 | sklearn.cluster |
| DBSCAN | 具有噪声的基于密度的聚类方法 | sklearn.cluster |
| BIRCH | 综合的层次聚类算法,可以处理大规模数据聚类 | sklearn.cluster |
这些不同的模型的使用方法大同小异,基本上是先用对应的函数建立模型,然后用.fit()方法训练模型,训练好之后可以用.label_方法给出样本数据的标签或者使用.predict()方法预测新的输入的标签.
下面介绍聚类结果可视化的工具-TSNE,定位是高维数据的可视化.一般输入的特征数是高维的(大于三维),难以直接以原特征对聚类结果进行展示,而TSNE提供一种有效的数据降维方式,让我们在二维或者3维空间展示聚类结果,我们使用TSNE对上述KMeans聚类结果以二维的方式展现出来:
#-*- coding: utf-8 -*-
#接k_means.py
from sklearn.manifold import TSNE
tsne = TSNE()
tsne.fit_transform(data_zs) #进行数据降维
tsne = pd.DataFrame(tsne.embedding_, index = data_zs.index) #转换数据格式
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号
#不同类别用不同颜色和样式绘图
d = tsne[r[u'聚类类别'] == 0]
plt.plot(d[0], d[1], 'r.')
d = tsne[r[u'聚类类别'] == 1]
plt.plot(d[0], d[1], 'go')
d = tsne[r[u'聚类类别'] == 2]
plt.plot(d[0], d[1], 'b*')
plt.show()
结果如图:
关联规则
关联规则也称为购物篮分析,最早为了发现购物车销售数据库中不同商品之间的关联关系,例如顾客可能购买了面包的顾客可能会购买牛奶,这就导致了一条关联规则"面包=>牛奶",其中面包称为规则前项,而牛奶称为后项,通过对面包降价售价进行促销,而适当提高牛奶售价,关联销售出的牛奶可能增加超市的整体利润.关联规则是数据挖掘中最活跃的研究方法之一,目的是在一个数据集中找出各项之间的关联关系,而这种关系并没有在数据中直接表示出来.
常用关联规则算法
| 算法名称 | 算法描述 |
|---|---|
| Apriori | 关联规则最常用也是最经典的挖掘频繁项集的算法,其核心是通过连接产生候选项及其支持度然后通过剪枝生成频繁项集 |
| FP-Tree | 针对Apriori算法的固有的多次扫描事务数据集的缺陷,提出的不产生候选频繁项集的方法,Apriori和FP-Tree都是寻找频繁项集的算法 |
| Eclat算法 | Eclat算法是一种深度优先算法,采用垂直数据表现形式,在概念格理论的基础上利用前缀的等价关系确定将搜索空间划分为较小的子空间 |
| 灰色关联法 | 分析和确定各因素之间的影响程度或是若干个子因素(子序列)对主因素(母序列)的贡献度而进行的一种分析方法 |
Apriori算法
以超市销售为例,提取关联规则最大的困难在于存在很多种商品时,可能商品的组合(规则的前项和后项)的数目会达到一种令人望而却步的程度.因而提取关联规则的算法需要从不同的方面入手,以减少可能的搜索空间的大小以及减少扫描的次数.Apriori是最经典的挖掘数据频繁集的算法,第一次在大数据集上可行的关联规则的提取.
- 关联规则和频繁项集
1.关联规则的一般形式
项集A,B同时发生的概率称为关联规则的支持度(也称相对支持度).
项集A发生,则项集B发生的概率为关联规则的置信度.
2.最小支持度和最小置信度
最小支持度是用户或专家定义的衡量支持度的一个阈值,表示项目集在统计意义上的最低重要性;最小置信度是用户或专家定义衡量置信度的一个阈值,表示关联规则最低可靠性,同时满足最小支持度阈值和最小置信度阈值的规则称为强规则.
3.项集
项集是项的集合.包含k个项的项集称为k项集,如集合(牛奶,麦片,糖)是一个三项集,项集的出现概率是包含所有项集的事务计数,又称为绝对支持度或支持度计数,如果项集I的相对支持度满足预定义的最小支持度阈值,则I是频繁项集,频繁项集通常记作k.
4.支持度计数
项集A的支持度计数是全部事务数据集中的包含项集A事务的个数的概率,称为项集的频率或者计数.已知项集的支持度计数,则规则A=>B的支持度和置信度很容易从所有事务计数,项集A和项集AUB的支持度计数推出.
举一个简单的关联规则例子,假设I是项的集合,即所有商品的无重复集合,给定一个交易数据库D,其中每个事务(Transaction)t是I的非空子集,即每一次购买都是所有商品的一部分,每一个交易都与一个唯一的标识符TID(Transaction ID)对应,关联规则在D中的支持度(support)是D中事务同时包含X、Y的百分比,即概率,置信度(confidence)是D中事物已经包含X的情况下,包含Y的百分比,即条件概率,如果满足最小支持度阈值和最小置信度阈值,这些阈值是根据挖掘需要人为设定,为频繁项集.
| TID | 网球拍 | 网 球 | 运动鞋 | 羽毛球 |
|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 0 |
| 2 | 1 | 1 | 0 | 0 |
| 3 | 1 | 0 | 0 | 0 |
| 4 | 1 | 0 | 1 | 0 |
| 5 | 0 | 1 | 1 | 1 |
| 6 | 1 | 1 | 0 | 0 |
用一个简单的例子说明,表1是顾客购买记录的数据库D,包含6个事务.项集I={网球拍,网球,运动鞋,羽毛球},考虑关联规则（频繁二项集）:网球拍与网球,事务1,2,3,4,6包含网球拍,事务1,2,6同时包含网球拍和网球,XY=3,D=6,支持度(XY)/D=0.5;X=5,置信度(X^Y)/X=0.6.若给定最小支持度α = 0.5,最小置信度β = 0.6,认为购买网球拍和购买网球之间存在关联.
Apriori算法代码为:
#apriori.py
#-*- coding: utf-8 -*-
from __future__ import print_function
import pandas as pd
#自定义连接函数，用于实现L_{k-1}到C_k的连接
def connect_string(x, ms):
x = list(map(lambda i:sorted(i.split(ms)), x))
l = len(x[0])
r = []
for i in range(len(x)):
for j in range(i,len(x)):
if x[i][:l-1] == x[j][:l-1] and x[i][l-1] != x[j][l-1]:
r.append(x[i][:l-1]+sorted([x[j][l-1],x[i][l-1]]))
return r
#寻找关联规则的函数
def find_rule(d, support, confidence, ms = u'--'):
result = pd.DataFrame(index=['support', 'confidence']) #定义输出结果
support_series = 1.0*d.sum()/len(d) #支持度序列
column = list(support_series[support_series > support].index) #初步根据支持度筛选
k = 0
while len(column) > 1:
k = k+1
print(u'\n正在进行第%s次搜索...' %k)
column = connect_string(column, ms)
print(u'数目：%s...' %len(column))
sf = lambda i: d[i].prod(axis=1, numeric_only = True) #新一批支持度的计算函数
#创建连接数据，这一步耗时、耗内存最严重。当数据集较大时，可以考虑并行运算优化。
d_2 = pd.DataFrame(list(map(sf,column)), index = [ms.join(i) for i in column]).T
support_series_2 = 1.0*d_2[[ms.join(i) for i in column]].sum()/len(d) #计算连接后的支持度
column = list(support_series_2[support_series_2 > support].index) #新一轮支持度筛选
support_series = support_series.append(support_series_2)
column2 = []
for i in column: #遍历可能的推理，如{A,B,C}究竟是A+B-->C还是B+C-->A还是C+A-->B？
i = i.split(ms)
for j in range(len(i)):
column2.append(i[:j]+i[j+1:]+i[j:j+1])
cofidence_series = pd.Series(index=[ms.join(i) for i in column2]) #定义置信度序列
for i in column2: #计算置信度序列
cofidence_series[ms.join(i)] = support_series[ms.join(sorted(i))]/support_series[ms.join(i[:len(i)-1])]
for i in cofidence_series[cofidence_series > confidence].index: #置信度筛选
result[i] = 0.0
result[i]['confidence'] = cofidence_series[i]
result[i]['support'] = support_series[ms.join(sorted(i.split(ms)))]
result = result.T.sort_values(['confidence','support'], ascending = False) #结果整理，输出
print(u'\n结果为：')
print(result)
return result
调用代码为:
#-*- coding: utf-8 -*-
#使用Apriori算法挖掘菜品订单关联规则
from __future__ import print_function
import pandas as pd
from apriori import * #导入自行编写的apriori函数
inputfile = '../data/menu_orders.xls'
outputfile = '../tmp/apriori_rules.xls' #结果文件
data = pd.read_excel(inputfile, header = None)
print(u'\n转换原始数据至0-1矩阵...')
ct = lambda x : pd.Series(1, index = x[pd.notnull(x)]) #转换0-1矩阵的过渡函数
b = map(ct, data.as_matrix()) #用map方式执行
data = pd.DataFrame(list(b)).fillna(0) #实现矩阵转换，空值用0填充
print(u'\n转换完毕。')
del b #删除中间变量b，节省内存
support = 0.2 #最小支持度
confidence = 0.5 #最小置信度
ms = '---' #连接符，默认'--'，用来区分不同元素，如A--B。需要保证原始表格中不含有该字符
find_rule(data, support, confidence, ms).to_excel(outputfile) #保存结果
输出结果:
结果为：
support confidence
e---a 0.3 1.000000
e---c 0.3 1.000000
c---e---a 0.3 1.000000
a---e---c 0.3 1.000000
c---a 0.5 0.714286
a---c 0.5 0.714286
a---b 0.5 0.714286
c---b 0.5 0.714286
b---a 0.5 0.625000
b---c 0.5 0.625000
a---c---e 0.3 0.600000
b---c---a 0.3 0.600000
a---c---b 0.3 0.600000
a---b---c 0.3 0.600000
其中e—a表示e发生推算a发生置信度为100%,支持度为30%,c—e---a表示ce同时发生能够推出a发生,置信度为100%,支持度为0.3.
- Apriori算法:使用候选集产生频繁项集
Apriori算法的主要思想是找出存在事务数据中最大频繁项集,在利用得到的最大频繁项集与预先设定的最小置信度阈值生成强关联规则.
1.Aprior的性质
频繁项集的非空子集也必须是频繁项集,根据该性质可以得出:向不是频繁项集I的项集添加事务A,新的项集IUA一定也不是频繁项集.
2.Aprior算法实现的两个过程
1.找出所有的频繁项集,在这个过程中连接步和剪枝步相互融合,最终得到最大的频繁项集Lk.
连接步:连接步的目的是找到K项集,对给定的最小支持度阈值,分别对I项候选集C1,剔除小于该阈值的项集得到I项的频繁集L1;下一步由L1自身连接产生的2项候选集C2,保留C2中满足约束条件的项集得到的到2项频繁集,记为L2,再下一步由L2与L3连接产生的3项候选集C3,保留C3中满足约束条件的项集得到3项频繁集,记为L3…这样循环下去,得到最大的频繁项集Lk.
若有两个k-1项集,每个项集按照"属性-值"(一般按值)的字母顺序进行排序.如果两个k-1项集的前k-2个项相同,而最后一个项不同,则证明它们是可连接的,即这个k-1项集可以联姻,即可连接生成k项集.使如有两个3项集:｛a, b, c｝{a, b, d},这两个3项集就是可连接的,它们可以连接生成4项集｛a, b, c, d｝.又如两个3项集｛a, b, c｝｛a, d, e｝,这两个3项集显示是不能连接生成3项集的.
剪枝步:剪枝步紧接着连接步,在产生候选项Ck的过程减少搜索空间的目的,若一个项集的子集不是频繁项集,则该项集肯定也不是频繁项集.很好理解,举一个例子,若存在3项集｛a, b, c｝,如果它的2项子集｛a, b｝的支持度即同时出现的次数达不到阈值,则｛a, b, c｝同时出现的次数显然也是达不到阈值的.因此,若存在一个项集的子集不是频繁项集,那么该项集就应该被无情的舍弃.
2.由频繁项集产生的强关联规则:由1已知未达到指定最小支持度阈值的项集已被剔除,如果剩下的这些规则又满足了最小置信度阈值,那么就挖掘出了强关联规则.
下面以餐饮行业数据为例,讲解Aprior关联规则算法挖掘实现过程,数据库部分点餐数据如图:
| 序列 | 时间 | 订单号 | 菜品ID | 菜品名称 |
|---|---|---|---|---|
| 1 | 2019/09/08 | 101 | 18491 | 街头肠粉 |
| … | … | … | … | … |
将表中的事务数据整理成关联规则模型所需的数据结构,从中抽取10个点餐订单作为事务数据集,设支持度为0.2(支持度计数为2),为方便起见,将菜品ID{18491,8842,8693,7794,8705} 分别简记为{a,b,c,d,e},见表
| 订单号 | 菜品ID | 订单号 | 菜品ID |
|---|---|---|---|
| 1 | a,c,e | 6 | b,c |
| 2 | b,d | 7 | a,b |
| 3 | b,c | 8 | a,b,c,e |
| 4 | a,b,c,d | 9 | a,b,c |
| 5 | a,b | 10 | a,c,e |
算法过程如图所示:
算法过程为:
过程一:找最大K项频繁集
1.算法简单的扫描所有的事务,事务中的每一项都是候选1项集的集合C1的成员,计算每一项的支持度,例如,p({a}) = 项集{a}的支持度计数/所有事务个数 = 7/10 =0.7
2.对C1各项的支持度与预先设定的最小支持度阈值进行比较,保留大于或等于该阈值的项,得1项频繁集L1.
3.扫描所有事务,L1与L1连接得候选2项集C2,并计算每一项的支持度.如P({a,b}) = 项集{a,b}的支持度计数/所有事务个数 = 5/10=0.5,接下来是剪枝,由于C2的每一项子集即L1都是频繁集,所以没有项集从C2中删除.
4.对C2中各项集的支持度与预先设定的最小支持度阈值进行比较,保留大于或等于该项阈值的项,得2项频繁集L2.
5.扫描所有事务,L2与L1连接得候选3项集C3,并计算每一项支持度,如P({a,b,c})=项集{a,b,c}的支持度计数/所有事务个数 = 3/10 =0.3,接下来的是剪枝步,L2与L1连接的所有的项集为:{a,b,c},{a,b,d},{a,b,e},{a,c,d},{a,c,e},{b,c,d},{b,c,e},根据算法,频繁集的所有非空子集也必须是频繁集,因为L2删去部分数据,所以{b,d},{b,e},{c,d}不包含在频繁子集L2中,即不是频繁集,应剔除,最后C3的项集只有{a,b,c},{a,c,e};
6.对C3的各项集的支持度与预先设定的最小支持度阈值进行比较,保留大于或等于的项,得3项频繁集L3.
7.L3项与L1项连接得候选4项集,剪枝后为空集,最后得最大3项频繁集为{a,b,c}和{a,c,e}.
由以上过程可知L1,L2,L3都是频繁集,最大项频繁集为L3.
过程二:由频繁集产生关联规则
置信度计算公式如下,计算A发生条件下B发生概率:
与python程序得到结果显然是一致的,就解释一条输出结果为例:客户同时点菜品A和B的概率为50%,点了菜品A,再点菜品B的概率是71.4%,知道了这些,就可以对客户进行智能推荐了,增加销量同时满足客户需求.
Rule Support Confidence
a---b 0.5 0.714286
时序模式
由于餐饮企业是生产和销售同步进行的,因此销售预测对餐饮企业是十分重要的,餐饮企业的销售预测可以看做基于时间序列的短期数据预测,预测对象为具体菜品销售量.常用按时间顺序排列的一组随机变量X1,X2,…,Xt来表示一个随机事件的时间序列,简记为{Xi};用x1,x2,…,xn或{xt,t = 1,2,…,n}表示该随机序列的n个有序观测值,称之为序列长度w为n的观测值序列.应用时间序列分析的目的就是给定一个已经被观测了的时间序列,预测该序列的未来值.
时间序列算法
常见的时间序列模型见图:
时间序列的处理
拿到一个观测值序列后,首先要对它的纯随机性和平稳性进行检验,这两个重要的检验称为序列的预处理,根据检验的结果可以将序列划分为不同的类型,对不同的类型的序列会采取不同的分析方法.
对于纯随机序列,又称白噪声序列,序列的各项之间没有任何相关关系,序列在进行着完全无序的随机波动,可以终止对该序列的分析,白噪声序列是没有任何信息可提取的平稳序列.对于平稳非白噪声序列,它的均值和方差是常数,现已有一套非常成熟的平稳序列的建模方法,通常是建立一个线性模型来拟合该序列的发展,借此提取该序列的有用信息.ARMA模型是最常用的平稳序列的拟合模型.对于非平稳序列,由于它的均值和方差不稳定,处理方法一般是将其转变为平稳序列,这样可以应用平稳时间序列的分析方法.如建立ARMA模型来进行相应的研究,如果一个时间序列经差分运算后具有平稳性,则该序列为差分平稳序列,可以使用ARIMA模型进行分析.
平稳性检验
平稳时间序列的定义
对于随机变量X,可以计算器均值(数学期望)μ,方差 σ 2 σ^2 σ2;对于两个随机变量X和Y,可以计算X,Y的协方差 c o v ( X , Y ) = E [ ( X − μ x ) ( Y − μ y ) ] cov(X,Y)=E[(X-μx)(Y-μy)] cov(X,Y)=E[(X−μx)(Y−μy)]和相关系数 ρ ( X , Y ) = c o v ( X , Y ) σ x σ y \rho (X,Y)=\frac{cov(X,Y)}{\sigma_x \sigma _y} ρ(X,Y)=σxσycov(X,Y),它们度量了两个不同事件之间的相互影响程度.
对于时间序列 X t , t ∈ T { X_{t},t\in T } Xt,t∈T,任意时刻的序列值 X t X_t Xt都是一个随机变量,每一个随机变量都会有均值和方差,记 X t X_t Xt的均值为 μ t \mu _t μt,方差为 σ t \sigma _{t} σt;任意取 t , s ∈ T t,s \in T t,s∈T,定义序列 X t { X_t } Xt的自协方差函数 γ ( t , s ) = E [ ( X t − μ t ) ( X s − μ s ) ] \gamma ( t,s )=E [ ( X_t-\mu _t )( X_s-\mu _s ) ] γ(t,s)=E[(Xt−μt)(Xs−μs)]和自相关系数 ρ ( t , s ) = c o v ( X t , X s ) σ t σ s \rho (t,s)=\frac{cov(X_t,X_s)}{\sigma _t\sigma _s} ρ(t,s)=σtσscov(Xt,Xs)(特别地, γ ( t , t ) = γ ( 0 ) = 1 , ρ 0 = 1 \gamma (t,t)=\gamma(0)=1,\rho _0=1 γ(t,t)=γ(0)=1,ρ0=1),之所以称它们为自协方差函数和自相关系数,是因为它们衡量的是同一个事件在两个不同时期(时刻t和s)之间的相关程度,就是度量自己过去的行为对自己现在的影响.
如果时间序列 X t , t ∈ T { X_{t},t\in T } Xt,t∈T在某一常数附近波动且波动范围有限,即有常数均值和方差,并且延迟k期的序列变量的自协方差和自相关系数是相等的或者说是延迟k期序列变量之间的影响程度是一样的,则称 X t , t ∈ T { X_{t},t\in T } Xt,t∈T为平稳序列.
平稳性的检验
对序列的平稳性检验有两种检验方法,一种是时序图和自相关图的特征作出判断的图检验,该方法操作简单,缺点是带有主观性,另一种是构造检验统计量进行检验的方法,目前最常用的是单位根检验.
时序图检验:根据平稳时间序列的均值和方差都为常数的性质,平稳序列的时序图显示该序列的值始终在一个常数附近随机波动,而且波动范围有界,如果有明显的周期性或者趋势性,那他通常不是平稳序列.
自相关图检验:平稳序列具有短期相关性,这个性质表明对平稳序列而言通常只有近期的序列值对现时值的影响比较明显,间隔越远的过去值对现时值影响最小,随着延迟期数k的增加,平稳序列的自相关系数 ρ k \rho _k ρk(延迟k期)会比较快的衰减趋向于零,并在零附近随机波动,而非平稳序列的自相关系数衰减的速度比较慢.
单位根检验:单位根检验是指检验序列中是否存在单位根,如果存在单位根就是非平稳时间序列.
纯随机性检验
如果一个序列是纯随机序列,那么它的序列值之间没有任何关系,即满足\gamma (k) = 0,k不为0这是一种理论上才会出现的理想状态,实际上纯随机序列的样本的自相关系数不会为0而是很接近0,并在0附近随机波动.
纯随机性检验也称白噪声检验,一般是构造检验统计量来检验序列的纯随机性,常用的检验统计量有Q统计量,LB统计量,由样本各延迟期数可以计算得到检验统计量,然后计算出对应p值,如果p值显著大于显著水平a,则表示该序列不能拒绝纯随机的原假设,可以停止对该序列的分析.
平稳时间序列分析
ARMA模型的全称是自回归移动平均模型,目前最常用的拟合平稳序列的模型,它又可细分为AR模型,MA模型和ARMA三大类,都可以看做多元线性回归模型.
AR模型
具有以下结构的模型称为p阶自回归模型,简记为AR§: X t = φ 0 + φ 1 X t − 1 + φ 2 X t − 2 + ⋯ + φ p X t − p + ε t X_{t}=\varphi_{0}+\varphi_{1}X_{t-1}+\varphi_{2}X_{t-2}+\cdots +\varphi_{p}X_{t-p}+\varepsilon _{t} Xt=φ0+φ1Xt−1+φ2Xt−2+⋯+φpXt−p+εt.
即在t时刻的随机变量 X t X_{t} Xt的取值是前p期 X t − 1 , X t − 2 , ⋯   , X t − p X_{t-1},X_{t-2},\cdots,X_{t-p} Xt−1,Xt−2,⋯,Xt−p的多元线性回归,认为 X t X_{t} Xt主要受过去p期的序列值的影响.误差项是当前的随机干扰 ε t \varepsilon _{t} εt,为零均值的白噪声序列.
平稳AR模型的性质见表:
| 统计量 | 性质 | 统计量 | 性质 |
|---|---|---|---|
| 均值 | 常数均值 | 自相关系数(ACF) | 拖尾 |
| 方差 | 常数方差 | 偏自相关系数(PACF) | p阶截尾 |
均值
对满足平稳性条件的AR§模型的方程,两边取期望,得:
E
(
X
i
)
=
E
(
φ
0
+
φ
1
X
t
−
1
+
φ
2
X
t
−
2
+
⋯
+
φ
p
X
t
−
p
+
ε
t
)
E(X_{i})=E(\varphi_{0}+\varphi_{1}X_{t-1}+\varphi_{2}X_{t-2}+\cdots +\varphi_{p}X_{t-p}+\varepsilon_{t})
E(Xi)=E(φ0+φ1Xt−1+φ2Xt−2+⋯+φpXt−p+εt),已知
E
(
X
i
)
=
μ
,
E
(
ε
t
)
=
0
E(X_{i})=\mu ,E(\varepsilon_{t})=0
E(Xi)=μ,E(εt)=0,所以解得
μ
=
φ
0
1
−
φ
1
−
φ
2
−
⋯
−
φ
p
\mu = \frac{\varphi_{0}}{1-\varphi_{1}-\varphi_{2}-\cdots -\varphi_{p}}
μ=1−φ1−φ2−⋯−φpφ0.
方差
平稳AR§方差有界,等于常数.
自相关系数(ACF)
平稳AR§模型的自相关系数
p
k
=
p
(
t
,
t
−
k
)
=
c
o
v
(
X
t
,
X
t
−
k
)
σ
t
σ
t
−
k
p_{k}=p(t,t-k)=\frac{cov(X_{t},X_{t-k})}{\sigma_t \sigma_{t-k}}
pk=p(t,t−k)=σtσt−kcov(Xt,Xt−k)呈指数的速度衰减,始终有非零取值,这个性质就是平稳AR§模型的自相关系数具有自相关性.
偏自相关系数(PACF)
对于一个平稳AR§模型求出延迟k期自相关系数
p
k
p_{k}
pk时,并未得到
X
t
X_{t}
Xt与
X
t
−
k
X_{t-k}
Xt−k之间单纯线性关系,
X
t
X_{t}
Xt同时受到中间k-1个变量影响,为了单纯测度
X
t
X_{t}
Xt对
X
t
−
k
X_{t-k}
Xt−k的影响,引入偏自相关系数.
平稳AR§模型的偏自相关系数具有p阶截尾性,这个性质与自相关系数拖尾性是AR§的重要识别依据.
MA模型
具有以下结构的模型称为q阶自回归模型,简记为MA(q): X t = μ + ε t − θ 1 ε t − 1 − θ 2 ε t − 2 − ⋯ − θ q ε t − q X_{t}=\mu + \varepsilon_{t}-\theta_{1}\varepsilon_{t-1}-\theta_{2}\varepsilon_{t-2}-\cdots-\theta_{q}\varepsilon_{t-q} Xt=μ+εt−θ1εt−1−θ2εt−2−⋯−θqεt−q,即在t时刻的随机变量 X t X_{t} Xt的取值 x t x_{t} xt是前q期的随机扰动 ε t − 1 , ε t − 2 , ⋯   , ε t − q \varepsilon_{t-1},\varepsilon_{t-2},\cdots,\varepsilon_{t-q} εt−1,εt−2,⋯,εt−q的多元线性函数,误差项是当期随机干扰 ε t \varepsilon_{t} εt,为零均值白噪声序列, μ \mu μ是序列 X t X_{t} Xt的均值,认为 X t X_{t} Xt主要是受过去q期的误差项的影响.平稳MA(q)模型性质见表:
| 统计量 | 性质 | 统计量 | 性质 |
|---|---|---|---|
| 均值 | 常数均值 | 自相关系数(ACF) | q阶截尾 |
| 方差 | 常数方差 | 偏自相关系数(PACF) | 拖尾 |
ARMA模型
具有以下结构的模型称为自回归移动平均模型,简记为ARMA(p,q):
x
t
=
ϕ
0
+
ϕ
1
x
t
−
1
+
ϕ
2
x
t
−
2
+
⋯
+
ϕ
q
x
t
−
p
+
ε
t
−
θ
1
ε
t
−
1
−
θ
2
ε
t
−
2
−
−
⋯
−
θ
q
ε
t
−
q
x_{t}=\phi_{0}+\phi_{1}x_{t-1}+\phi_{2}x_{t-2}+\cdots+\phi_{q}x_{t-p}+\varepsilon_{t}-\theta_{1}\varepsilon_{t-1}-\theta_{2}\varepsilon_{t-2}--\cdots-\theta_{q}\varepsilon_{t-q}
xt=ϕ0+ϕ1xt−1+ϕ2xt−2+⋯+ϕqxt−p+εt−θ1εt−1−θ2εt−2−−⋯−θqεt−q,
即在t时刻的随机变量
X
t
X_{t}
Xt的取值
x
t
x_{t}
xt是前p期
x
t
−
1
,
x
t
−
2
,
⋯
 
,
x
t
−
p
x_{t-1},x_{t-2},\cdots,x_{t-p}
xt−1,xt−2,⋯,xt−p和前q期
ε
t
−
1
,
ε
t
−
2
,
⋯
 
,
ε
t
−
q
\varepsilon_{t-1},\varepsilon_{t-2},\cdots,\varepsilon_{t-q}
εt−1,εt−2,⋯,εt−q的多元线性函数,误差项是当前随机干扰
ε
t
\varepsilon_{t}
εt,为零均值白噪声序列,认为
x
t
x_{t}
xt主要是受过去p期的序列值和过去q期的误差项的共同影响.
特别的,当q=0时,是AR§模型;当p=0时,是MA(q)模型.平稳ARMA模型性质见表:
| 统计量 | 性质 | 统计量 | 性质 |
|---|---|---|---|
| 均值 | 常数均值 | 自相关系数(ACF) | 拖尾 |
| 方差 | 常数方差 | 偏自相关系数(PACF) | 拖尾 |
平稳时间序列建模
某个时间序列经过预处理,被判定为平稳非白噪声序列,就可以使用ARMA模型进行建模.计算出平稳非白噪声序列
X
t
{X_{t}}
Xt的自相关系数和偏自相关系数,再由AP§模型,MA(q)模型和ARMA(p,q)的自相关系数和偏自相关系数性质选择合理模型,平稳时间序列建模步骤如图:
计算ACF和PACF,之后进行ARMA模型识别,也称,模型订阶,由AR§模型,MA(q)模型和ARMA(p,q)的自相关系数和偏自相关系数,选择合适模型,识别原则见表:
| 模型 | 自相关系数(ACF) | 偏自相关系数(PACF) |
|---|---|---|
| AR§ | 拖尾 | p阶截尾 |
| MA(q) | q阶截尾 | 拖尾 |
| ARMA(p,q) | p阶拖尾 | q阶拖尾 |
下一步估计模型中未知参数的值并进行参数进行检验,模型检验,模型优化,模型应用(进行短期预测).
非平稳时间序列分析
对非平稳时间序列分析方法可以分为确定性因素分解的时序分析和随机分析两大类.确定性因素分析把所有序列的变化归为四个因素(长期因素,季节变动,循环变动和随机波动)的综合影响,其中长期趋势和季节变动的规律信息通常比较容易提取,而随机因素导致的波动则非常难以确定和分析,对随机信息浪费严重,会导致模型拟合精度不够理想.
随机时序分析根据时间序列的不同特点,随机时序分析可以建立的模型有ARIMA模型,残差自回归模型,季节模型,异方差模型等.
差分运算
p阶差分:相距一期的两个序列值之间的减法运算称为1阶差分运算.
k步差分:相距k期的两个序列值之间的减法运算称为k步差分运算.
ARIMA模型
差分运算具有强大的确定性信息提取能力,许多非平稳序列差分后会显示出平稳序列的性质,这时称这个非平稳序列为差分平稳序列,对差分平稳序列可以使用ARMA模型进行拟合,ARIMA模型的实质就是差分运算与ARMA模型的组合,差分平稳时间序列建模步骤如图所示;
下面对demo/data/arima_data.xls中2015/1/1-2015/2/16某餐厅的销售数据进行建模:
#-*- coding: utf-8 -*-
#arima时序模型
import pandas as pd
#参数初始化
discfile = '../data/arima_data.xls'
forecastnum = 5
#读取数据，指定日期列为指标，Pandas自动将“日期”列识别为Datetime格式
data = pd.read_excel(discfile, index_col = u'日期')
#时序图
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号
data.plot()
plt.show()
#自相关图
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data).show()
#平稳性检测
from statsmodels.tsa.stattools import adfuller as ADF
print(u'原始序列的ADF检验结果为：', ADF(data[u'销量']))
#返回值依次为adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore
#差分后的结果
D_data = data.diff().dropna()
D_data.columns = [u'销量差分']
D_data.plot() #时序图
plt.show()
plot_acf(D_data).show() #自相关图
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(D_data).show() #偏自相关图
print(u'差分序列的ADF检验结果为：', ADF(D_data[u'销量差分'])) #平稳性检测
#白噪声检验
from statsmodels.stats.diagnostic import acorr_ljungbox
print(u'差分序列的白噪声检验结果为：', acorr_ljungbox(D_data, lags=1)) #返回统计量和p值
from statsmodels.tsa.arima_model import ARIMA
data[u'销量'] = data[u'销量'].astype(float)
#定阶
pmax = int(len(D_data)/10) #一般阶数不超过length/10
qmax = int(len(D_data)/10) #一般阶数不超过length/10
bic_matrix = [] #bic矩阵
for p in range(pmax+1):
tmp = []
for q in range(qmax+1):
try: #存在部分报错，所以用try来跳过报错。
tmp.append(ARIMA(data, (p,1,q)).fit().bic)
except:
tmp.append(None)
bic_matrix.append(tmp)
bic_matrix = pd.DataFrame(bic_matrix) #从中可以找出最小值
p,q = bic_matrix.stack().idxmin() #先用stack展平，然后用idxmin找出最小值位置。
print(u'BIC最小的p值和q值为：%s、%s' %(p,q))
model = ARIMA(data, (p,1,q)).fit() #建立ARIMA(0, 1, 1)模型
model.summary2() #给出一份模型报告
model.forecast(5) #作为期5天的预测，返回预测结果、标准误差、置信区间。
检验序列的平稳性
获得原始序列时序图(日期与销量)
原始序列的自相关图:
获得原始序列的单位根检验:
原始序列的ADF检验结果为： (1.813771015094526, 0.9983759421514264, 10, 26, {'1%':
-3.7112123008648155, '5%': -2.981246804733728, '10%': -2.6300945562130176}, 299
.4698986602418)
从时序图我们可以看出该序列有明显的单调递增趋势,可以判断为非平稳序列,自相关系数并没有趋向0,说明序列有长期相关性,单位根统计量0.99验证表明对应p显著大于0.05,最终将该序列判断为非平稳序列(非平稳序列一定不是白噪声序列).
对原始序列进行一阶差分,并进行平稳性和白噪声检验
一阶差分之后的序列的时序图:
一阶差分之后的自相关图:
一阶差分之后的单位根检验:
差分序列的ADF检验结果为： (-3.1560562366723537, 0.022673435440048798, 0, 35, {'1
%': -3.6327426647230316, '5%': -2.9485102040816327, '10%': -2.6130173469387756},
287.5909090780334)
结果显示,一阶差分之后的序列的时序图在均值附近波动,自相关图具有很强的短期相关性,单位根检验小于0.05,所以一阶差分之后的序列是平稳序列.
对一阶差分之后的序列进行白噪声检验:
差分序列的白噪声检验结果为： (array([11.30402222]), array([0.00077339]))
输出的p值远小于0.05,所以一阶差分之后的序列是平稳非白噪声序列.
对一阶差分之后的平稳非白噪声序列拟合ARMA模型
下面进行模型定阶,模型定阶就是确定p和q.
第一种就是人为识别的方法,一阶差分后自相关图显示一阶截尾,偏自相关图显示拖尾,所以考虑用MA模型拟合1阶差分后的序列,即对原始序列建立ARIMA(0,1,1)模型.
第二种就是相对最优模型识别,计算ARMA(p,q).当p和q值均小于等于3的所有组合的BIC信息量,取其中BIC信息量达到最小的模型阶数,计算完成BIC矩阵如下:
经过计算可得BIC最小值是p,q的值:
BIC最小的p值和q值为：0、1
所以对原始序列建立ARIMA(0,1,1)模型,虽然两种方法建立模型是一样的,但是模型并非唯一,可以检验ARIMA(1,1,0)和ARIMA(1,1,1),这两个模型也可以通过检验.
对一阶差分后的序列拟合AR(1)模型进行分析,可以得到一份模型参数:
Coef. Std.Err. t P>|t| [0.025 0.975]
----------------------------------------------------------------------
const 49.9564 20.1390 2.4806 0.0182 10.4847 89.4281
ma.L1.D.销量 0.6710 0.1648 4.0712 0.0003 0.3480 0.9941
-----------------------------------------------------------------------------
ARIMA模型预测
应用ARIMA(0,1,1)对表中销售数据做为期5天的预测,预测结果为:
(array([4873.9667493 , 4923.92317644, 4973.87960359, 5023.83603073,
5073.79245787]), array([ 73.08574293, 142.32679918, 187.542821 , 223.80
81869,
254.95704265]), array([[4730.72132537, 5017.21217324],
[4644.96777602, 5202.87857687],
[4606.30242887, 5341.4567783 ],
[4585.19056646, 5462.48149499],
[4574.08583666, 5573.49907907]]))
需要说明的是预测时期越长,预测误差越大,这是时间预测的典型特点.
Python主要时序模式算法
Python实现时序模式主要库是StatsModels,算法模型主要是ARIMA模型,在使用该模型时,需要进行一系列判别操作,主要包括平稳性检验,白噪声检验,是否差分,AIC和BIC指标值,模型定阶,最后再做预测,其主要相关函数见表:
| 函数名 | 函数功能 | 所属工具箱 |
|---|---|---|
| acf() | 计算相关系数 | statsmodels.tsa.stattools |
| plot_acf() | 画自相关系数图 | statsmodels.graphics.tsaplots |
| pacf() | 计算偏相关系数 | statsmodels.tsa.stattools |
| plot_pacf() | 画偏相关系数图 | statsmodels.graphics.tsaplots |
| adfuller() | 对观测值序列进行单位根检验 | statsmodels.tsa.stattools |
| diff() | 对观测值序列进行差分运算 | Pandas对象自带的方法 |
| ARIMA() | 创建一个ARIMA时序模型 | statsmodels.tsa.arima_model |
| summary()或summaty2 | 给出一份ARIMA模型的报告 | ARIMA模型对象自带的方法 |
| aic/bic/hqic | 计算ARIMA模型的AIC/BIC/HQIC指标值 | ARIMA模型对象自带的变量 |
| forecast() | 应用构建的时序模型进行预测 | ARIMA模型对象自带的方法 |
| acorr_ljungbox() | Ljung-Box检验,检验是否为白噪声 | statsmodels.stats.diagnostic |
acf()
功能:计算自相关系数
使用格式:autocorr=acf(data,unbiased=False,nlags=40,qstat=False,fft=False,alpha=None)
输入参数data为观测值序列(即为时间序列,可以是DataFrame或Series),返回参数autocorr为观测值序列自相关函数,其余为可选参数,如qstat=True,返回Q统计量和对应p值.
plot_acf()
功能:画自相关系数图
使用格式:p = plot_acf(data)
返回一个Matplotlib对象,可以使用.show()方法显示图像.
pacf()/plot_pacf()
功能:计算偏相关系数/画偏相关系数图
使用格式:使用跟acf()/plot_acf()类似.
adfuller()
功能:对观测值序列进行单位根检验(ADF test)
使用格式:h = adfuller(Series,maxlag=None,regression=‘c’,autolag=‘AIC’,store=False,regresults=False)
输入参数Series为一维观测序列值,返回值依次为adf,pvalue,usedlag,nobs,critical values,icbest,regresults,resstore.
diff()
功能:对观测值序列进行差分运算
使用格式:D.diff(),D为Pandas的DataFrame或Series.
arima
功能:设置时序模式的建模参数,创建ARIMA时序模型
使用格式:arima = ARIMA(data,(p,l,q)).fit()
data参数为输入的时间序列,p,q为对应的阶,d为差分次数.
summary()/summary2()
功能:生成已有模型的报告
使用格式:arima.summary()或arima.summary2()
其中,arima为已经建立好的ARIMA模型,返回一份格式化的模型报告,包含模型的系数,标准误差,p值,AIC和BIC等详细指标.
aic/bic/hqic
功能:计算ARIMA模型的AIC,BIC,HQIC指标值
使用格式:arima.aic/arima.bic/arima.hqic
其中,arima为已经已经建立好的ARIMA模型,返回的是Model时序模型得到的AIC,BIC,HQIC指标值.
forecast()
功能;用得到的时序模型进行预测
使用格式:a,b,c = arima.forecast(num)
输入参数num为要预测的天数,arima为已经建立好的ARIMA模型.a为返回的预测值,b为预测的误差,c为预测置信区间.
acorr_ljungbox()
功能:检测是否为白噪声序列
使用格式:acorr_ljungbox(data,lags=1),
输入参数data为时间序列数据,lags为滞后数,返回统计量和p值.
离群点检测
就餐饮企业而言,经常会碰到以下问题:如何根据客户的消费记录检测是否为异常刷卡消费,如何检测是否有异常订单?这两类问题可以通过离群点检测来解决.
离群点检测是数据挖掘中重要的一部分它的任务是发现与大部分其他对象显著不同的对象,大部分数据挖掘方法都将这种差异信息视为噪声而丢弃,然而在一些应用中,罕见的数据可能蕴含着更大的研究价值.在数据散布图中,离群点远离其他数据点,因为离群点的属性值明显偏离期望的或常见的属性值,所以离群点检测也称偏差检测.
离群点的主要成因有:数据来源于不同的类,自然变异,数据测量和收集误差.离群点的类型大致分类如下表:
| 分类标准 | 分类名称 | 分类描述 |
|---|---|---|
| 从数据范围 | 全局离群点和局部离群点 | 从整体来看,某些对象没有离群特征,但是从局部来看,却显示了一定的离群特性如图所示,C是全局离群点,D是局部离群点 |
| 从数据类型 | 数值型离群点和分类型离群点 | 这是以数据集的属性类型进行划分的 |
| 从属性个数 | 一维离群点和多维离群点 | 一个对象可能有一个或多个属性 |
离群点检测方法
常见的离群点检测方法如下表:
| 离群点检测方法 | 方法描述 | 方法评估 |
|---|---|---|
| 基于统计 | 构建一个概率分布模型,并计算对象符合该模型的概率,把具有低概率的对象视为离群点 | 基于统计模型的离群点检测方法的前提是必须知道数据集服从什么分布;对于高维数据检验效果比较差 |
| 基于邻近度 | 通常可以在数据对象之间定义邻近性度量,把远离大部分点的对象视为离群点 | 简单的二维三维数据可以选做散点图观查,大数据集不适用,对参数选择敏感,不能处理具有不同密度的数据集 |
| 基于密度 | 考虑数据集可能存在不同密度区域这一事实,从基于密度的观点分析,离群点是在低密度区域的对象,一个对象离群点得分是该对象周围密度的逆 | 给出了对象是离群点的定量度量,并且即使数据具有不同的区域也能够很好的处理;大数据集不适用 |
| 基于聚类 | 利用聚类检测离群点的方法是丢弃远离其他簇的小簇,另一种更系统的方法,首先聚类所有对象,然后评估对象属于簇的程度(离群点得分) | 聚类算法产生的簇的质量对该算法产生的离群点的质量影响非常大 |
基于模型的离群点检测方法
通过估计概率分布的参数来建立一个数据模型,如果一个数据对象不能很好的与该模型拟合,则它可能不服从该分布,它是一个离群点.
一元正态分布离群点检测
正态分布是统计学最常用的分布,若随机变量x的密度函数
φ
(
x
)
=
1
2
π
e
−
(
x
−
μ
)
2
2
σ
2
(
x
∈
R
)
\varphi(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}} (x\in R)
φ(x)=2π1e−2σ2(x−μ)2(x∈R),则称x服从正态分布
N
(
μ
,
σ
)
N(\mu,\sigma)
N(μ,σ),其中参数
μ
\mu
μ和
σ
\sigma
σ分布为均值和标准差.图显示N(0,1)的密度函数.
N(0,1)的数据对象出现在该分布两边的尾部的机会很小,因此可以用它作为检测数据对象是否是离群点的基础.数据对象落在3倍标准差中心区域之外的概率仅有0.0027.
混合模型的离群点检测
混合模型是一种特殊的统计模型,它使用若干统计分布对数据建模,每一个分布对应一个簇,而每个分布的参数提供对应簇的描述,通常使用中心和发散描述.混合模型的数据看做从不同的概率分布得到的观测值的集合.概率分布可以是任何分布,但是通常是多元正态的,因为这种类型的分布不难理解,容易从数学处理,并且已经证明在许多情况下都能产生良好的结果,类型分布可以使用椭圆簇建模.
混合模型数据的产生过程为:给定几个类型相同但参数不同的分布,随机的选取一个分布并由它产生一个对象,重复该过程m次,其中m是对象的个数.具体的讲,假定有K个分布和m个对象为
χ
=
{
x
1
,
x
2
,
⋯
 
,
x
m
}
\chi={\left\{x_{1},x_{2},\cdots ,x_{m}\right\}}
χ={x1,x2,⋯,xm}.设第j个分布的参数为
a
j
a_{j}
aj,并设A是所有参数的集合,即
A
=
{
a
1
,
a
2
,
⋯
 
,
a
k
}
A=\left\{a_{1},a_{2},\cdots,a_{k}\right\}
A={a1,a2,⋯,ak}.,则
P
(
x
i
∣
a
j
)
P(x_{i}|a_{j})
P(xi∣aj)是第i个对象来自第j个分布的参数的概率.选取第j个分布产生一个对象的概率由权值
ω
j
(
1
⩽
j
⩽
K
)
\omega_{j}(1\leqslant j\leqslant K)
ωj(1⩽j⩽K)给定,其中权值(概率)受限于其和为1的约束,即
∑
j
=
1
K
ω
j
=
1
\sum_{j=1}^{K}\omega_{j}=1
∑j=1Kωj=1.于是,对象x的概率由以下公式给出:
P
(
x
∣
A
)
=
∑
j
=
1
K
ω
j
P
j
(
x
∣
θ
j
)
P(x\mid A)=\sum_{j=1}^{K}\omega_{j}P_{j}(x|\theta_{j})
P(x∣A)=∑j=1KωjPj(x∣θj).如果对象以独立的方式产生,则整个对象集的概率是每个个体对象
x
j
x_{j}
xj的概率乘积,公式如下:
对于混合模型,每一个分布描述一个不同的组,即一个不同的簇.通过统计的方法,可以由数据估计这些分布的参数.也可以识别到那个对象属于哪个簇.然而,混合模型只是给出具体对象属于特定簇的概率,聚类时,混合模型方法假定数据来自混合概率分布并且每个簇可以用这些分布之一识别.同样的对于离群点检测,用两个分布的混合模型建模,一个分布为正常数据,一个为离群点.
聚类和离群点检测的目标都是估计分布的参数,以最大化数据的总似然:先将所有数据对象放入正常数据集,这是离群点为空集;再用一个迭代过程将数据对象从正常数据集转移到离群点集,该转移可以提高数据的总似然.具体操作如下:
假设数据集U包含来自两个概率分布的数据对象:M是大多数正常数据对象分布,而N是离群点数据对象分布,数据总概率为
U
(
x
)
=
(
1
−
λ
)
M
(
x
)
+
λ
N
(
x
)
U(x)=(1-\lambda)M(x)+\lambda N(x)
U(x)=(1−λ)M(x)+λN(x),其中x是一个数据对象,
λ
∈
[
0
,
1
]
\lambda\in [0,1]
λ∈[0,1],给出离群点的期望比例.分布M由数据估计得到,而分布N通常取均匀分布.设
M
t
M_{t}
Mt和
N
t
N_{t}
Nt分别为时刻t正常数据和离群点数据的集合,初始t=0,
M
0
M_{0}
M0=D,而
N
0
≠
ϕ
N_{0}\neq \phi
N0̸=ϕ.根据混合模型公式推导,在整个数据集的似然和对数似然可分别由下面两式给出:
其中
P
D
P_{D}
PD,
P
M
t
P_{M_{t}}
PMt,
P
N
t
P_{N_{t}}
PNt分别是D,
M
t
M_{t}
Mt和
N
t
N_{t}
Nt的概率分布函数.
因为正常数据集对象的数量比离群点对象数量大很多,因此当一个数据对象移动到离群点集后,正常数据对象分布变化不大,在这种情况下,此外假设离群点服从均匀分布,则移动到离群点集每一个数据对象对离群点的似然贡献一个固定的量,这样当一个数据对象移动到离群点集时,数据总似然的改变粗略的等于该数据对象在均匀分布下的概率(用
λ
\lambda
λ加权)减去该数据对象在正常数据点分布概率(用
1
−
λ
1-\lambda
1−λ加权).从而离群点由这样一些数据对象组成,这样数据对象在均匀分布下的概率比正常数据对象分布下概率要高.
基于聚类的离群点检测
聚类分析常用来发现局部强相关的对象组,而异常检测用来发现不与其他对象强相关的对象,本文主要介绍两种基于聚类的离群点检测方法.
丢弃远离其他簇的小簇
一种利用聚类检测离群点的方法是丢弃远离其他簇的小簇,通常,通过该过程可以简化为丢弃小于某个最小阈值的所有簇.
需要最小簇的大小和小簇与其他簇之间的距离的阈值,而且这种方案对簇个数的选择高度敏感,使用这个方案很难将离群点得分附加到对象上.在图中,聚类的簇数K=2,可以直接看到其中包含5个对象的小簇远离大部分对象,可以设为离群点.
基于原型的聚类
首先聚类所有对象,然后评估对象属于簇的程度(离群点得分),使用对象到它的簇中心的距离来度量属于簇的程度,特别的,如果删除一个对象使得该目标显著改进,则可以将该对象视为离群点.例如,在K均值算法中,删除原离其相关簇中心的对象能够显著改进该簇的误差平方和(SSE).
基于原型的聚类,主要有两种方法评估对象属于簇的程度(离群点得分):一是度量对象到簇原型的距离,并用它作为对象的离群点得分;二是考虑簇具有不同的密度,可以使用簇到原型的相对距离,相对距离是点到质心的距离与簇中所有点到质心距离的中位数之比.如下图,如果选择聚类簇数K=3,则对象ABC应分别属于距离它们最近的簇,但相对于簇内的其他对象,这3个点又分别远离各自的簇,所以有理由怀疑对象A,B,C是离群点.
诊断步骤如下:进行聚类.选择聚类算法(如K-Means算法),将样本集聚为K簇,并找到各簇的质心;计算各对象到它最近质心的距离;计算各对象到它的最近的质心的相对距离;与给定的阈值作比较.
对demo/data/consumption_data.xls数据进行聚类分析,计算各个样本到各自聚类中心的距离,分析离群样本,得到离散点检测距离误差图,代码:
#-*- coding: utf-8 -*-
#使用K-Means算法聚类消费行为特征数据
import numpy as np
import pandas as pd
#参数初始化
inputfile = '../data/consumption_data.xls' #销量及其他属性数据
k = 3 #聚类的类别
threshold = 2 #离散点阈值
iteration = 500 #聚类最大循环次数
data = pd.read_excel(inputfile, index_col = 'Id') #读取数据
data_zs = 1.0*(data - data.mean())/data.std() #数据标准化
from sklearn.cluster import KMeans
model = KMeans(n_clusters = k, n_jobs = 4, max_iter = iteration) #分为k类，并发数4
model.fit(data_zs) #开始聚类
#标准化数据及其类别
r = pd.concat([data_zs, pd.Series(model.labels_, index = data.index)], axis = 1) #每个样本对应的类别
r.columns = list(data.columns) + [u'聚类类别'] #重命名表头
norm = []
for i in range(k): #逐一处理
norm_tmp = r[['R', 'F', 'M']][r[u'聚类类别'] == i]-model.cluster_centers_[i]
norm_tmp = norm_tmp.apply(np.linalg.norm, axis = 1) #求出绝对距离
norm.append(norm_tmp/norm_tmp.median()) #求相对距离并添加
norm = pd.concat(norm) #合并
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号
norm[norm <= threshold].plot(style = 'go') #正常点
discrete_points = norm[norm > threshold] #离群点
discrete_points.plot(style = 'ro')
for i in range(len(discrete_points)): #离群点做标记
id = discrete_points.index[i]
n = discrete_points.iloc[i]
plt.annotate('(%s, %0.2f)'%(id, n), xy = (id, n), xytext = (id, n))
plt.xlabel(u'编号')
plt.ylabel(u'相对距离')
plt.show()
距离误差图为:
分析距离误差图可知,如果设置距离阈值为2,那么所给的数据中有8个离散点,在聚类时应将数据剔除.
扫一扫
1762
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
