看书标记——R语言
【R语言数据分析项目精解:理论、方法、实战】
Chapter 4 指标监控系统
4.1项目背景、目标方案
4.1.1项目背景
指标过多,人工跟踪工作量大,效率低,所以希望可以针对重要指标制定一套自动报警机制,把每天有异常的指标自动输出,达到有效降低人力成本的目的,进一步的,可以在这方面通过历史数据对未来进行预测。
4.1.2项目目标
针对重要指标建立预测模型,通过预测模型的95%预测上下限建立监控范围。这样就可以一举两得,预测模型可以对未来进行预估,另外95%上下限建立的范围可以用于监控,若当天数值超出当天预测值的监控范围,则报警。最后,通过可视化工具前端展示整个需求就可以了。
4.1.3项目方案
分成两种情况,节假日和非节假日。非节假日建立平稳时期指标量化模型,节假日建立影响数据量化模型,然后加一个示性函数,得到包含节假日的综合量化模型。
4.2项目技术理论简介
4.2.1时间序列基本统计量
时间序列均值和方差、自协方差函数和自相关函数
4.2.2数据观测与描述性统计
1.折线图:时间序列数据最直观的统计图形,能很好地反映整个趋势的走向和周期性,从而为下一步的分析提供指引性的先验知识。
2.平稳性:(严平稳)当序列所有的统计性质都不会随着时间的推移而发生变化;(宽平稳)当序列的统计性质只与时间间隔的变化而发生变化,以上两种都说明序列稳定。
3.检验方法:正态性检验(序列服从正态分布,一定平稳,所以也叫平稳性检验,此时H0:数据服从正态分布);自相关图(自相关系数是否随着延迟期数的增加而迅速衰减到零,通常平稳序列具有6期以内的短期相关性,而后平稳);单位根检验(单位根检验有单位根说明序列非平稳,非平稳序列在d阶差分后成平稳序列,称序列d阶单整。DF检验用于一阶自相关模型检验,ADF检验用于高阶自相关、高阶移动平均的模型,PP检验用于扰动项之间存在序列相关,KPSS检验用于序列是平稳或趋势平稳,优于DF检验)
4.非平稳序列平稳化:不通过平稳性检验的就变换序列使得变化后的序列通过平稳性检验。常用的变换方法有:差分变换、对数变换、Box-Cox变换。
4.2.3随机性
序列平稳后需对数据的随机性进行检验,若数据是随机的,那么对数据的近一步探讨建模将会是没有意义的。
Barlett定理、LB随机检验统计量
4.2.4周期性
除平稳性和随机性外,周期性也是时间序列的一个基本属性,除了折线图外,周期性还可以在谱分析的周期图中查看,高峰值可计算周期。
4.2.5节假日模式识别
针对带有周期的时间序列数据,提出周期相同天的配对t检验识别算法,找出带有周期的时间序列数据节假日影响模式的识别方法,具体步骤如下:
1.取原始数据》2.标注可疑的异常日期(主观标记)》3.周期相同天分组》4.周期相同天分组均值和》5.配对样本t检验
4.2.6建模数据集
1.挑选建模数据集
根据行业特征选取适合进行数据分析的时间范围内 的数据。
2.寻找离群点
噪声数据会对指标序列产生不利影响,一般来说,指标的不稳定性主要有:框架程序出错、导致记录不全;ETL调度工具运行不稳,导致用户行为记录不全甚至缺失;业务本身导致指标异常(举办活动等);节假日因素。
先对“每周相同天定义”,同周期的第一天作为一个数据集,对每个集合使用检验方法(Grubbs检验、Dixon检验(小样本)、线箱图法(没有分布约束)和基于密度的局部异常因子法),挑出离群点。
3.平滑离群点
对于被剔除的离群点,通常会选择均值进行填补,具体的填补公式见P123详解。
4.2.7指标监控方法(不含节假日)
1.ARIMA模型
- AR模型:自回归模型又称为中心化的AR§模型,平稳过程为因果平稳( AR模型的自相关系数是呈复指数衰减– 有拖尾性;AR模型的偏自相关系数有截尾性)
- MA模型:移动平均模型,MA模型总是平稳的,只有MA模型是依历史可逆的,该过程才具有预测意义(自相关系数q阶截尾;偏自相关系数q阶拖尾)
- ARMA模型:AR和MA的结合,所以性质与其相似
- 差分运算:p阶差分(p次只跨一个时刻);k步差分(一次跨k个时刻)
- ARIMA模型:应用差分运算提取确定性信息后建模,所以ARIMA就是ARMA与差分运算的结合。
ARIMA模型模式识别
ACF法和PACF法判断p、q的值,或者AIC、BIC信息准则选取模型(autoarima中常用)
残差分析(真实值与模型预测值比较)
标准残差——时间图:零线上下无规则波动
标准残差的独立性检验:ACF图里滞后阶相关系数都在标准差范围内,则认为满足独立性检验。
残差随机性检验:广义方差检验和LB统计量检验(残差具有随机性则表明模型将规律部分全部提取)
正态性检验:Shapiro-Wilk检验和QQ图判断数据是否正态(残差符合正态分布,模型效果较好)
2.季节乘积ARIMA模型
季节乘积ARMA模型和季节乘积ARIMA模型(是否有差分运算的区别)
3.指数平滑
最近的数据对未来的影响更大,故应赋予较大的权重。指数平滑法本质上是对过去的数据进行加权和,距今越近,权重越高。
一次指数平滑:适用于序列无明显趋势,缺点是序列出现线性趋势时,一次指数平滑会出现严重滞后。
二次指数平滑:适用于序列存在线性趋势,缺点是序列出现周期性,则二次平滑预测序列偏差较大。
三次指数平滑:适用于序列存在线性及周期性趋势。累加和累乘,当季节变化在整个系列中大致恒定时,优选加法,而当季节变化与系列水平成比例变化时,优选乘法。
指数平滑对历史数据采取加权的方式比较符合实际情况,模型容易搭建,且具有自适应性,会根据数据变化自动变换。但是指数平滑模型对时间序列的转折点洞悉不足,适合短期预测。
4.质量控制图
用于监控产品零件生产品质
(1)适用条件:当监控指标变化较平稳且近似正态时,也可以用于互联网指标监控中。
(2)均值-极差控制图(hat_x—MR图)
(3)过程能力指数:衡量了过程加工内在的一致性,表现的是加工过程固有的波动状态。
5.判断异常值的方法
季节乘积ARIMA模型和指数平滑模型由样本内拟合均方误差和样本外预测均方误差SSE的值选择,若拟合和预测选择的不一致,则可以,计算样本内两拟合均方误差的倒数得到两模型权重,随后以权重整合模型,得到预测值和置信限。
4.2.9R语言实例代码
1.Grubbs检验
######################################################################
#函数功能:Grubbs检验
#参数说明:x:要进行判断的数据
#####################################################################
library(outliers)
grubbs<-function(x){
x<-round(x,4)
grubbs_outliers<-c()
grubbs_p.value<-c()
grubbs_g.value<-c()
grubbs_g<-c()
grubbs_minormax<-c()
grubbs_pvalue<-c()
#设定初始值判断值
grubbs_p<-0
while(grubbs_p<0.05){
grubbs_outliers<-c(grubbs_outliers,grubbs_minormax)
grubbs_p.value<-c(grubbs_p.value,grubbs_pvalue)
grubbs_g<-c(grubbs_g,grubbs_g.value)
#去除最大、最小值然后再进行检验
if(sum(x==grubbs_minormax)!=0)x<-x[-which(x==grubbs_minormax)]
#若数据各值都相等则跳出
if(sd(x)==0) break
#进行grubbs检验
grubbs_test<-grubbs.test(x,type=10,opposite=F,two.sided=F)
grubbs_p<-grubbs_test$p.value
grubbs_pvalue<-grubbs_test$p.value
grubbs_g.value<-grubbs_test$statistic[1]
grubbs_a<-strsplit(grubbs_test$alternative," ",fixed=T)
grubbs_minormax<-as.numeric(unlist(grubbs_a)[3])
}
outliner_res<-data.frame(outliers=grubbs_outliers,gvalue=grubbs_g,pvalue=grubbs_p.value)
return(outliner_res)
}
tt<-c(8.3,5.5,14,7.5,4.7,9,6.5,10.2,7.7,6.2)
grubbs(tt)
代码说明:grubbs. test是Grubbs检验实现函数,其调用格式为Grubbs. test (x,type=10,opposite=FALSE,two. sided=FALSE)。其中,参数x为检验数据向量;type=10表示检验一个异常值,type=11表示检验分别处于两个端点的异常值,type=20表示检验两个在同一侧的异常值;two. sided表示双边检验。上述代码通过Grubbs. test函数找出一个异常值,然后去除这个异常值再次检验,直到没有异常值为止,最后输出所有检验出来的异常值、G统计量和pvalue。
2.Dixon检验
######################################################################
#函数功能:Dixon检验
#参数说明:x:要进行判断的数据
#####################################################################
library(outliers)
dixon<-function(x){
dixon_p.value<-c()
dixon_q.value<-c()
dixon_q<-c()
dixon_pvalue<-c()
dixon_outliers<-c()
dixon_minormax<-c()
dixon_p<-0
while(dixon_p<0.05){
dixon_outliers<-c(dixon_outliers,dixon_minormax)
dixon_p.value<-c(dixon_p.value,dixon_pvalue)
dixon_q<-c(dixon_q,dixon_q.value)
if(sum(x==dixon_minormax)!=0)x<-x[-which(x==dixon_minormax)]
if(sd(x)==0) break
dixon_test<-dixon.test(x,type=0,opposite=F,two.sided=F)
dixon_p<-dixon_test$p.value
dixon_pvalue<-dixon_test$p.value
dixon_q.value<-dixon_test$statistic[1]
dixon_a<-strsplit(dixon_test$alternative," ",fixed=T)
dixon_minormax<-as.numeric(unlist(dixon_a)[3])
}
outliner_res<-data.frame(outliers=dixon_outliers,qvalue=dixon_q,pvalue=dixon_p.value)
return(outliner_res)
}
tt<-c(8.3,5.5,14,7.5,4.7,9,6.5,10.2,7.7,6.2)
dixon(tt)
代码说明:函数思路和Grubbs检验的函数思路相似。dixon. test为Dixon检验的主要函数,其调用格式为dixon. test (x,type=10,
参数x为检验数据;
opposite=FALSE,two. sided=TRUE)。其中,
type=10表示检验数据集为8~10个的数据,type=21表示检验数据集为11~13个的数据,type=22表示检验数据集为14个或14个以上的数据;two. sided表示双边检验。通过自定义的Dixon检验函数,没有发现异常值。
3.LOF检验
########################################################################
#函数功能:lof检验
#参数说明:x:要进行判断的数据
########################################################################
library(DMwR)
library(outliers)
loffun<-function(klb,kup,indexnum){
lof_value<-c()
indexdata<-indexnum
indexname<-"indexdata"
index_lof<-data.frame(indexdata=indexdata)
for (k in klb:kup){
lof_score<-paste("lof_score",k,"<-lofactor(",indexname,",",k,")",sep="")
add_lof<-paste("index_lof$lof_score",k,"<-lof_score",k,sep="")
eval(parse(text=lof_score))
eval(parse(text=add_lof))
cr_lofresult<-paste("lof_value<-apply(cbind(lof_value,lof_score",k,"),1,max)",sep="")
eval(parse(text=cr_lofresult))
}
index_lof$lofvalue<-round(lof_value,2)
index_lof$num<-1
lof_thres<-grubbs(index_lof$lofvalue)
index_lof$lofflag<-c(rep(0,nrow(index_lof)))
index_lof[index_lof$lofvalue%in%lof_thres$outliers,]$lofflag<-1
lof_result<-subset(index_lof,lofflag==1)
lof_result[order(-lof_result$lofvalue),]
lof_data<-lof_result[,c("indexdata","lofvalue")]
lof_data[order(-lof_data$lofvalue),]
return(lof_data)
}
klb<-10
kup<-14
#生成测试数据
indexdata<-c(2,1,2,3,2,3,4,5,6,5,4,3,22,1,2,3,4,0,2,33)
loffun(klb,klb,indexdata)
代码说明: lofactor函数用于LOF离群点的判断,其中k为每个族中最少包含的元素个数。通过循环函数,对不同k的可能取值依次计算LOF值,然后取最大的一个作为该数值的LOF值;接着用前面的Grubbs函数确定LOF阈值(找出离群点),最终LOF值离群点对应的数值即为样本离群点。假定k取10~14,通过循环检验,最终找出了10个离群点。
4.指数平滑
HoltWinters(co2) #三次指数平滑
5.均值-极差图
#####################################################################
#函数功能:极差均值控制图
#参数说明:x:要进行判断的数据
#####################################################################
xmr<-function(index){
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,1))
color<-c("green","red")
title1<-paste("cr_xmr_单值",sep="")
title2<-paste("cr_xmr_移动极差",sep="")
#计算极差均值数据
R<-abs(diff(index))
Rbar<-mean(R)
xbar<-mean(index)
xucl<-xbar+2.66*Rbar
xlcl<-xbar-2.66*Rbar
Rucl<-3.267*Rbar
Rlcl<-0
#作图
plot(index,ylim=c(xlcl-0.2*xlcl,xucl+0.2*xucl),type='o',pch=19,main=title1,col=color[1])
abline(h=xucl,col=color[2])
abline(h=xlcl,col=color[2])
plot(R,ylim=c(Rlcl-0.2*Rlcl,Rucl+0.2*Rucl),type='o',pch=19,main=title2,col=color[1])
abline(h=Rucl,col=color[2])
abline(h=Rlcl,col=color[2])
par(opar)
interval<-list(ucl=xucl,lcl=xlcl)
return(interval)
}
#生成指标
indexdata<-c(0.1071,0.1097,0.1069,0.116,0.0893,0.0877,0.1181,0.1107,0.1122)
xmr(indexdata)
4.3项目实践
4.3.1数据概览
文中用到的是订单支付转化率,CR是一个相对平稳的指标,结合了订单和流量这两个核心指标,也是KPI重要的核心指标。
1.折线图
#数据获取
setwd("C:\\Users\\用户路径")
crtot<-read.csv("2011~2016.10总体数据.csv",header=TRUE)
crtot$orderdate<-as.Date(crtot$orderdate)
indexdata<-crtot$indexdata
x.text<-crtot$orderdate
#总体折线图
plot(indexdata,type="l",xaxt="n",xlab="日期",ylab="数值")
abline(h=mean(indexdata))
axis(1,at=1:length(indexdata),labels=x.text,tick=FALSE)
#每年折线图
crtot$year<-as.numeric(format(crtot$orderdate,"%Y"))
flt_11<-subset(crtot,year==2011)
flt_12<-subset(crtot,year==2012)
flt_13<-subset(crtot,year==2013)
flt_14<-subset(crtot,year==2014)
flt_15<-subset(crtot,year==2015)
flt_16<-subset(crtot,year==2016)
id=c("2011年","2012年","2013年","2014年","2015年","2016年")
col=c("black","red","orange","purple","yellow","blue")
plot(flt_11$indexdata,type='l',ylim=c(0,0.2),col=col[1],main="各年数据")
legend("topleft",legend=id,horiz=T,pch=15,col=col,cex=0.8,bty="n")
lines(flt_12$indexdata,col=col[2])
lines(flt_13$indexdata,col=col[3])
lines(flt_14$indexdata,col=col[4])
lines(flt_15$indexdata,col=col[5])
lines(flt_16$indexdata,col=col[6])
每年同期的折线图基本平稳
2.平稳性检验
library(TSA)
adfresult<-adf.test(indexdata,alt="stationary") #adf检验
adfresult
3.随机性检验
Box.test(indexdata, type="Ljung-Box",lag=6)
Box.test(indexdata, type="Ljung-Box",lag=12)
4.周期性检验
根据行业情况判断周期,普遍为7天,文中有用谱分析中的周期图来验证7天的假设是否正确。每年的指标数据大致相同;节假日对数据整体影响较大,所以先提出固定节假日的数据。
flt_pfx<-subset(flt_16,orderdate>='2016-06-01' & orderdate<='2016-09-09') #选取样本数据
pr<-periodogram(flt_pfx$indexdata) #画出谱分析中的周期图
score1<-pr$freq[which(pr$spec==max(pr$spec))] #计算最高点对应的周期
f<-1/score1 ;f #周期
由最高尖峰的频率来计算,1/0.1388=7.2046天
4.3.2节假日模式识别
1.节假日模式分类
根据各指标节前、节中、节后的不同变化趋势,对节假日期间的指标进行大致分类。
2.转化率指标节假日模式识别
文中用的是国庆节相关数据,首先要对指标数据打上假日标识,大致确定节前、节中、节后的时间段。
library(plyr)
holidayflag<-read.csv("holiday.csv",header=TRUE,stringsAsFactors=FALSE) #读取节假日标识数据
holidayflag$orderdate<-as.Date(holidayflag$orderdate)
flt<-read.csv("2011~2016.10总体数据.csv",header=TRUE) #读取指标数据
flt$orderdate<-as.Date(flt$orderdate)
flt<-join(flt,holidayflag) #将节假日标识打在指标数据上
flt$year<-as.numeric(format(flt$orderdate,"%Y"))
#每年数据分别存储
flt_11<-subset(flt,year==2011)
flt_12<-subset(flt,year==2012)
flt_13<-subset(flt,year==2013)
flt_14<-subset(flt,year==2014)
flt_15<-subset(flt,year==2015)
flt_16<-subset(flt,year==2016)
#国庆前后明细数据图
id=c("2011年","2012年","2013年","2014年","2015年","2016年")
col=c("black","red","orange","purple","yellow","blue")
xtext<-paste(format(flt_11_2$orderdate,"%m"),format(flt_11_2$orderdate,"%d"),sep='-')
plot(flt_11_2$indexdata,type='l',xaxt='n',col=col[1],ylim=c(0,0.2),main="国庆前后明细数据")
axis(1,at=1:61,labels=xtext,tick=FALSE)
legend("topleft",legend=id,horiz=T,pch=15,col=col,cex=0.8,bty="n")
lines(flt_12_2$indexdata,col=col[2])
lines(flt_13_2$indexdata,col=col[3])
lines(flt_14_2$indexdata,col=col[4])
lines(flt_15_2$indexdata,col=col[5])
lines(flt_16_2$indexdata,col=col[6])
然后周期相同天组成各自对应的样本集,计算出均值与之匹配,实现样本的配对
##周期相同天数据整理
flt_11$number<-c(rep(1:7,nrow(flt_11)%/%7),c(1:7)[1:(nrow(flt_11)%%7)])
flt_12$number<-c(rep(1:7,nrow(flt_12)%/%7),c(1:7)[1:(nrow(flt_12)%%7)])
flt_13$number<-c(rep(1:7,nrow(flt_13)%/%7),c(1:7)[1:(nrow(flt_13)%%7)])
flt_14$number<-c(rep(1:7,nrow(flt_14)%/%7),c(1:7)[1:(nrow(flt_14)%%7)])
flt_15$number<-c(rep(1:7,nrow(flt_15)%/%7),c(1:7)[1:(nrow(flt_15)%%7)])
flt_16$number<-c(rep(1:7,nrow(flt_16)%/%7),c(1:7)[0:(nrow(flt_16)%%7)])
#----------离群点判断----------------
library(outliers)
library(DMwR)
##Grubbs'Test
Grubbs<-function(initial,cycle){
grubbs<-function(x){
x<-round(x,6)
grubbs_outliers<-c()
grubbs_minormax<-c()
grubbs_p<-0
while(grubbs_p<0.05){
grubbs_outliers<-c(grubbs_outliers,grubbs_minormax)
if(sum(x==grubbs_minormax)!=0)x<-x[-which(x==grubbs_minormax)]
if(sd(x)==0) break
grubbs_test<-grubbs.test(x,type=10,opposite=F,two.sided=F)
grubbs_p<-grubbs_test$p.value
grubbs_a<-strsplit(grubbs_test$alternative," ",fixed=T)
grubbs_minormax<-as.numeric(unlist(grubbs_a)[3])
}
return(grubbs_outliers)
}
Grubbs.outliers<-tapply(initial$indexdata,initial$number,grubbs)
initial$Grubbs.test<-rep(0,nrow(initial))
for(i in 1:cycle){
initial$Grubbs.test[which(initial$indexdata%in%Grubbs.outliers[[i]]&initial$number==i)]<-1
}
return(initial)
}
##Boxplot
Boxplot<-function(initial,cycle){
Boxplot.outliers<-tapply(initial$indexdata,initial$number,function(x) boxplot.stats(x)$out)
initial$Boxplot.test<-rep(0,nrow(initial))
for(i in 1:cycle){
initial$Boxplot.test[which(initial$indexdata%in%Boxplot.outliers[[i]]&initial$number==i)]<-1
}
return(initial)
}
##Lofactor
Lof<-function(initial){
initial$Lof.test<-rep(0,nrow(initial))
Lof.scores<-lofactor(initial$indexdata,k=floor(nrow(initial)/10))
Lof.outliers<-initial$indexdata[which(Lof.scores>1.75)]
initial$Lof.test[which(Lof.scores>1.75)]<-1
return(initial)
}
holvalue<-function(tbdata){
res_after4<-subset(tbdata,orderdate>=paste(unique(tbdata$year),"-06-01",sep='') & orderdate<=paste(unique(tbdata$year),"-09-16",sep=''))
grubbs<-Grubbs(res_after4,7)
boxplot<-Boxplot(res_after4,7)
lof<-Lof(res_after4)
res_after4$Grubbs.test<-grubbs$Grubbs.test
res_after4$Boxplot.test<-boxplot$Boxplot.test
res_after4$Lof.test<-lof$Lof.test
res_after4$abnflag<-res_after4$Grubbs.test+res_after4$Boxplot.test+res_after4$Lof.test
tb_adj<-subset(res_after4,abnflag<2)
tb_adjmean<-aggregate(tb_adj$indexdata,by=list(tb_adj$number),mean)
names(tb_adjmean)[1]<-"number"
tb_hol<-subset(tbdata,complete.cases(tbdata$flag))
tb_rb<-merge(tb_hol,tb_adjmean,all.tb_hol=TRUE)
return(tb_rb)
}
############相同周期点配对数据输出调整格式#############
flt_11hol_adj<-flt_11hol[order(flt_11hol$flag),]
flt_12hol_adj<-flt_12hol[order(flt_12hol$flag),]
flt_13hol_adj<-flt_13hol[order(flt_13hol$flag),]
flt_14hol_adj<-flt_14hol[order(flt_14hol$flag),]
flt_15hol_adj<-flt_15hol[order(flt_15hol$flag),]
flt_16hol_adj<-flt_16hol[order(flt_16hol$flag),]
flt_11hol_adj$date<-paste(format(flt_11hol_adj$orderdate,'%m'),'-',format(flt_11hol_adj$orderdate,'%d'),sep='')
flt_12hol_adj$date<-paste(format(flt_12hol_adj$orderdate,'%m'),'-',format(flt_12hol_adj$orderdate,'%d'),sep='')
flt_13hol_adj$date<-paste(format(flt_13hol_adj$orderdate,'%m'),'-',format(flt_13hol_adj$orderdate,'%d'),sep='')
flt_14hol_adj$date<-paste(format(flt_14hol_adj$orderdate,'%m'),'-',format(flt_14hol_adj$orderdate,'%d'),sep='')
flt_15hol_adj$date<-paste(format(flt_15hol_adj$orderdate,'%m'),'-',format(flt_15hol_adj$orderdate,'%d'),sep='')
flt_16hol_adj$date<-paste(format(flt_16hol_adj$orderdate,'%m'),'-',format(flt_16hol_adj$orderdate,'%d'),sep='')
names(flt_11hol_adj)[c(3,6)]<-c("2011index","2011avg")
names(flt_12hol_adj)[c(3,6)]<-c("2012index","2012avg")
names(flt_13hol_adj)[c(3,6)]<-c("2013index","2013avg")
names(flt_14hol_adj)[c(3,6)]<-c("2014index","2014avg")
names(flt_15hol_adj)[c(3,6)]<-c("2015index","2015avg")
names(flt_16hol_adj)[c(3,6)]<-c("2016index","2016avg")
flt_11hol_adj<-flt_11hol_adj[,c("flag","date","2011index","2011avg")]
flt_12hol_adj<-flt_12hol_adj[,c("flag","date","2012index","2012avg")]
flt_13hol_adj<-flt_13hol_adj[,c("flag","date","2013index","2013avg")]
flt_14hol_adj<-flt_14hol_adj[,c("flag","date","2014index","2014avg")]
flt_15hol_adj<-flt_15hol_adj[,c("flag","date","2015index","2015avg")]
flt_16hol_adj<-flt_16hol_adj[,c("flag","date","2016index","2016avg")]
flt_tot_adj<-merge(flt_11hol_adj,flt_12hol_adj)
flt_tot_adj<-merge(flt_tot_adj,flt_13hol_adj)
flt_tot_adj<-merge(flt_tot_adj,flt_14hol_adj)
flt_tot_adj<-merge(flt_tot_adj,flt_15hol_adj)
flt_tot_adj<-merge(flt_tot_adj,flt_16hol_adj)
write.csv(flt_tot_adj,"相同周期点配对数据.csv") #结果输出到本地
先用两种及以上方法剔除对均值会有影响的离群点,然后计算均值后的配对样本数据调整格式并输出。
#------3)配对T检验
flt_tot<-rbind(flt_11hol,flt_12hol,flt_13hol,flt_14hol,flt_15hol,flt_16hol)
#配对T检验函数
tvalue<-function(id,flt_tot){
aa<-subset(flt_tot,flag==id)$indexdata
bb<-subset(flt_tot,flag==id)$x
t_value<-t.test(aa,bb,paired=TRUE)$p.value
return(t_value)
}
#针对每个日期,进行配对T检验
ttest_pvalue<-unlist(lapply(sort(unique(flt_tot$flag)),tvalue,flt_tot))
result_pvalue<-data.frame(flag=sort(unique(flt_12hol$flag)),pvalue=ttest_pvalue)
#计算均值
indexdata_avg<-aggregate(flt_tot$indexdata,by=list(flt_tot$flag),mean)
x_avg<-aggregate(flt_tot$x,by=list(flt_tot$flag),mean)
names(x_avg)[2]<-"xavg"
x_avg_df<-merge(indexdata_avg,x_avg,all.indexdata_avg=TRUE)
names(x_avg_df)[1]<-'flag'
ttest_result<-merge(x_avg_df,result_pvalue)
ttest_result$x_cf<-rep("higher",nrow(ttest_result))
ttest_result[which(ttest_result$x<=ttest_result$xavg),]$x_cf<-"lower"
ttest_result$pvalueflag<-rep("F",nrow(ttest_result))
ttest_result[which(ttest_result$pvalue<0.01),]$pvalueflag<-"T"
4.3.3模型数据集的建立
1.挑选数据集
cr<-read.csv("建模数据1031-带异常值.csv",header=TRUE) #读取数据
cr$orderdate<-as.Date(cr$orderdate)
cr<-subset(cr,cr$orderdate<='2016-09-11')
cycle<-7 #设置周期
cr$number<-rep(1:cycle,len=nrow(cr))
xtext<-paste(format(cr$orderdate,"%m"),format(cr$orderdate,"%d"),sep='-') #画出折线图
plot(cr$indexdata,type='l',xaxt='n')
axis(1,at=1:length(cr$orderdate),labels=xtext,tick=FALSE)
2.寻找离群点
#--Grubbs检验
Grubbs.outliers<-tapply(cr$indexdata,cr$number,grubbs)
cr$grubbsflag<-c(rep(0,nrow(cr)))
for(i in 1:cycle){
if(nrow(Grubbs.outliers[[i]])>0){
grubbs_df<-subset(Grubbs.outliers[[i]],pvalue<=0.1)
grubbs_outliers<-grubbs_df$outliers
cr$grubbsflag[which(cr$indexdata%in%grubbs_outliers&cr$number==i)]<-1
}
}
head(cr)
#--Dixon检验
Dixon.outliers<-tapply(cr$indexdata,cr$number,dixon)
cr$dixonflag<-c(rep(0,nrow(cr)))
for(i in 1:cycle){
if(nrow(Dixon.outliers[[i]])>0){
dixon_df<-subset(Dixon.outliers[[i]],pvalue<=0.1)
dixon_outliers<-dixon_df$outliers
cr$dixonflag[which(cr$indexdata%in%dixon_outliers&cr$number==i)]<-1
}
}
head(cr)
#--线箱图
Boxplot.outliers<-tapply(cr$indexdata,cr$number,function(x) boxplot.stats(x)$out)
cr$boxplotflag<-rep(0,nrow(cr))
for(i in 1:cycle){
cr$boxplotflag[which(cr$indexdata%in%Boxplot.outliers[[i]]&cr$number==i)]<-1
}
head(cr)
#--LOF局部离群点
klb<-10
kup<-14
crlof<-loffun(klb,kup,cr$indexdata)
cr$lofflag<-c(rep(0,nrow(cr)))
cr$lofflag[which(cr$indexdata%in%crlof$indexdata)]<-1
head(cr)
#--整合四种方法
cr$outlierflag<-cr$grubbsflag+cr$dixonflag+cr$boxplotflag+cr$lofflag
#查看总体结果
head(cr[order(-cr$outlierflag),],10)
至少有两种方法判定为离群点(“1”)的数据可以是真正的离群点。
3.离群点平滑:普遍采用均值填补
4.3.4指标监控(非节假日)
1.数据概览
建模
crdf<-read.csv("建模数据1031-带异常值 平滑后 adj.csv",header=TRUE) #获取数据
crdf$orderdate<-as.Date(crdf$orderdate)
#9月11日之前用来建模,9月12日~9月18日的用来检验
cr_df<-subset(crdf,orderdate<='2016-09-11')
cr_tot<-cr_df$indexdata
plot(cr_tot,type='o') #画折线图
abline(h=mean(cr_tot))
#计算自相关图和偏自相关
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2))
acf(cr_tot,lag.max=60)
pacf(cr_tot,lag.max=60)
par(opar)
#随机性检验
Box.test(cr_tot, type="Ljung-Box",lag=6)
Box.test(cr_tot, type="Ljung-Box",lag=12)
#平稳化检验
library(TSA)
adfresult<-adf.test(cr_tot,alt="stationary") #adf检验
原始序列具有7天的周期;Ljung-Box检验,p<0.05拒绝H0:序列随机的假设,有相关性信息;ADF检验,p<0.05拒绝H0:序列有单位根的假设,序列平稳。
2.乘积季节模型的建立
#选择ARIMA(0,1,2)×(2,0,1)7.
library(forecast)
arimatt_opt<-arima(cr_tot,order=c(0,1,2),seasonal=list(order=c(2,0,1),period=7),method='ML')
arimatt_opt
#残差分析
res<-arimatt_opt$res #计算残差
#残差-时间图
plot(rstandard(arimatt_opt),ylab='Standardized Residuals',main='标准残差-时间')
abline(h=0)
acf(res,main="ACF of residuals",lag.max=60) #残差独立性检验
#随机性检验
#----广义方差检验
library(portes)
plot(gvtest(arimatt_opt,1:60)[,4],main="广义方差检验",ylab='p-value ',xlab="lag",pch=16,ylim=c(0,1))
abline(h=0.05,lty=2)
#----LB检验
Box.test(res, type="Ljung-Box",lag=6)
#正态性检验
qqnorm(res)
qqline(res)
sol<-shapiro.test(res) #原假设:数据服从正太分布
sol
#---模型拟合与预测
prop.fore = predict(arimatt_opt, n.ahead =7)
#95%预测上线限
u<- prop.fore$pred + 1.96* prop.fore$se
l<-prop.fore$pred-1.96*prop.fore$se
cr_ts<-ts(cr_tot,frequency=1)
ts.plot(cr_ts, prop.fore$pred, col=c("black","blue"))
lines(u, col="red", lty="dashed")
lines(l, col="red", lty="dashed")
#样本内拟合
fit<-cr_tot-res
lines(fit,col="chocolate")
#样本外预测
cryz<-subset(crdf,orderdate>='2016-09-12' & orderdate<='2016-09-18')
realvalue<-cryz$indexdata
outlinepre_ar<-data.frame(pre=prop.fore$pred,realvalue=realvalue)
3.三次指数平滑
#模型拟合
cr_ts <- ts(cr_tot, frequency=7,start=c(1987,1))
plot.ts(cr_ts)
cr_ht <- HoltWinters(cr_ts)
plot(cr_ht)
#残差分析
e<-cr_ht$x-cr_ht$fitted[,1] #计算残差
plot(e,ylab='Standardized Residuals',main='标准残差-时间') #残差-时间图
abline(h=0)
acf(e,main="ACF of residuals",lag.max=60) #独立性检验
Box.test(e, type="Ljung-Box",lag=6) #随机性检验
Box.test(e, type="Ljung-Box",lag=12)
#正态性检验 #QQ图
qqnorm(e)
qqline(e)
sol<-shapiro.test(e) #原假设:数据服从正太分布
#(3)拟合与预测
library("forecast")
fitvalue<-fitted(cr_ht)[,1] #样本内预测
plot(cr_ts)
lines(fitvalue,col='chocolate')
cr_ht_fore <- forecast.HoltWinters(cr_ht, h=7) #样本外预测
plot.forecast(cr_ht_fore)
#95%样本外预测上下限
l<-cr_ht_fore[5]$lower[,2]
u<-cr_ht_fore[6]$upper[,2]
#样本外预测均方误
outlinepre_hw<-data.frame(realvalue=realvalue,prevalue=cr_ht_fore[[4]])
#(4)整合乘积季节模型和Holt-Winters模型
#乘积季节模型样本内拟合均方误
inline_mse_ar<-mean(abs(arimatt_opt$res)^2)
#Holt-Winter样本内拟合均方误
inline_mse_hw<-mean(e^2)
#乘积季节模型样本外拟合均方误
outline_mse_ar<-mean((outlinepre_ar$pre-outlinepre_ar$realvalue)^2)
#Holt-Winter样本外拟合均方误
outline_mse_hw<-mean((outlinepre_hw$realvalue-outlinepre_hw$prevalue)^2)
modelcf<-matrix(rep(0,4),nrow=2)
rownames(modelcf)<-c("样本内拟合均方误","样本外预测均方误")
colnames(modelcf)<-c("季节乘积模型","Holt-Winters")
modelcf[1,]<-c(inline_mse_ar,inline_mse_hw)
modelcf[2,]<-c(outline_mse_ar,outline_mse_hw)
4.质量监控图
#(1)对每周相同天集合做正态性检验
crdf<-read.csv("建模数据1031-带异常值 平滑后 adj.csv",header=TRUE)
crdf$orderdate<-as.Date(crdf$orderdate)
cycle<-7
cr_df<-subset(crdf,orderdate<='2016-09-11')
cr_df$number<-rep(1:cycle,len=nrow(cr_df))
#正态性检验函数
normtest<-tapply(cr_df$indexdata,list(cr_df$number),shapiro.test)
weekday<-c("周一","周二","周三","周四","周五","周六","周日")
normpvalue<-c()
for (i in 1:7){
normpvalue<-c(normpvalue,normtest[[i]]$p.value)
}
normtest<-data.frame(weekday=weekday,pvalue=normpvalue)
normtest
#(2)制作X-MR分析用控制图
##################################################################
#函数功能:极差均值控制图
#参数说明:df 数据框
#i:整数值
##################################################################
xmr<-function(df,i){
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,1))
color<-c("green","red")
index<-subset(df,number==i)$indexdata
title1<-paste("cr_xmr_单值:第",i,"组",sep="")
title2<-paste("cr_xmr_移动极差:第",i,"组",sep="")
R<-abs(diff(index))
Rbar<-mean(R)
xbar<-mean(index)
xucl<-xbar+2.66*Rbar
xlcl<-xbar-2.66*Rbar
Rucl<-3.267*Rbar
Rlcl<-0
plot(index,ylim=c(xlcl-0.2*xlcl,xucl+0.2*xucl),type='o',pch=19,main=title1,col=color[1])
abline(h=xucl,col=color[2])
abline(h=xlcl,col=color[2])
plot(R,ylim=c(Rlcl-0.2*Rlcl,Rucl+0.2*Rucl),type='o',pch=19,main=title2,col=color[1])
abline(h=Rucl,col=color[2])
abline(h=Rlcl,col=color[2])
par(opar)
interval<-list(ucl=xucl,lcl=xlcl)
return(interval)
}
xmr_1<-xmr(cr_df,1)
#均值和移动极差都在范围内,对有超出范围的部分,需要对该组数据进行分析和修正
xmr_3<-xmr(cr_df,3)
index<-subset(cr_df,number==3)$indexdata
abs(diff(index))
##所以需要去除第12个值后再进行计算
index<-subset(cr_df,number==3)$indexdata
abs(diff(index))
i<-3
index<-index[-length(index)]
title1<-paste("cr_xmr_单值:第",i,"组",sep="")
title2<-paste("cr_xmr_移动极差:第",i,"组",sep="")
R<-abs(diff(index))
Rbar<-mean(R)
xbar<-mean(index)
xucl_3_adj<-xbar+2.66*Rbar
xlcl_3_adj<-xbar-2.66*Rbar
Rucl<-3.267*Rbar
Rlcl<-0
color<-c("green","red")
plot(index,ylim=c(xlcl_3_adj-0.2*xlcl_3_adj,xucl_3_adj+0.2*xucl_3_adj),type='o',pch=19,main=title1,col=color[1])
abline(h=xucl_3_adj,col=color[2])
abline(h=xlcl_3_adj,col=color[2])
plot(R,ylim=c(Rlcl-0.2*Rlcl,Rucl+0.2*Rucl),type='o',pch=19,main=title2,col=color[1])
abline(h=Rucl,col=color[2])
abline(h=Rlcl,col=color[2])
##由质量监控图知,第三组均值和移动极差均在范围内
#(3)上下限确定
xmr_1<-xmr(cr_df,1)
xmr_2<-xmr(cr_df,2)
xmr_4<-xmr(cr_df,4)
xmr_5<-xmr(cr_df,5)
xmr_6<-xmr(cr_df,6)
xmr_7<-xmr(cr_df,7)
weekday<-c("周一","周二","周三","周四","周五","周六","周日")
ucl_interval<-c(xmr_1$ucl,xmr_2$ucl,xucl_3_adj,xmr_4$ucl,xmr_5$ucl,xmr_6$ucl,xmr_7$ucl)
lcl_interval<-c(xmr_1$lcl,xmr_2$lcl,xlcl_3_adj,xmr_4$lcl,xmr_5$lcl,xmr_6$lcl,xmr_7$lcl)
xmr_interval<-data.frame(weekday=weekday,ucl=round(ucl_interval,4),lcl=round(lcl_interval,4))
##给出转化率指标的最小值后,可以计算得到过程能力指数,将之与登记表对比可得过程稳定。
#### 4.3.5 节假日指标监控
##### 1.预测节假日数据
```r
library(TSA)
library(portes)
library(forecast)
crdf<-read.csv("建模数据1031-带异常值 平滑后 adj.csv",header=TRUE)
crdf$orderdate<-as.Date(crdf$orderdate)
cr_df<-subset(crdf,orderdate<='2016-09-18')
cr_tot<-cr_df$indexdata
arimatt_opt<-arima(cr_tot,order=c(0,1,2),seasonal=list(order=c(2,0,1),period=7),method='ML')
prop.fore <- predict(arimatt_opt, n.ahead =25)
#95%预测上线限
ucl_ar<- prop.fore$pred + 1.96* prop.fore$se
lcl_ar<-prop.fore$pred-1.96*prop.fore$se
cr_hol<-subset(crdf,orderdate>='2016-09-19' & orderdate<='2016-10-13')
crhol<-data.frame(indexdata=cr_hol$indexdata,pred=prop.fore$pred)
crhol$diff<-crhol$indexdata-crhol$pred
crhol
2.差值建模
cr_tot_diff<-crhol$diff
plot(cr_tot_diff,type='l')
pr<-periodogram(cr_tot_diff)
prspec<-sort(pr$spec,decreasing=TRUE)
f_1<-pr$freq[which(pr$spec==prspec[1])] ##第一个周期频率
f_2<-pr$freq[which(pr$spec==prspec[2])] ##第二个周期频率
3.模型拟合
#生成衍生变量,谱分析拟合
gydf<-data.frame(id=c(1:length(cr_tot_diff)),indexdata=cr_tot_diff)
gydf$v1_cos<-cos(2*pi*f_1*gydf$id)
gydf$v1_sin<-sin(2*pi*f_1*gydf$id)
gydf$v2_cos<-cos(2*pi*f_2*gydf$id)
gydf$v2_sin<-sin(2*pi*f_2*gydf$id)
#线性拟合
tt<-lm(formula=indexdata ~ v1_cos+v1_sin+v2_sin,data=gydf)
summary(tt) ##R^2接近1,拟合效果较好
#拟合图
plot(as.vector(cr_tot_diff),type='l')
lines(fitted(tt),col="blue")
4.残差分析
#残差图
res<-tt$residuals #计算残差
plot(tt,which=1) #残差图 近似水平,无规律
library(car)
durbinWatsonTest(tt) #DW独立性检验,p-value>0.05,符合独立性假设。
#同方差性
library(car)
ncvTest(tt) ##p-value>0.05,符合同方差性。
shapiro.test(res) #正态性检验
#预测上下限
#节假日第一天预测值
id_pre<-1
v1_cos_pre1<-cos(2*pi*f_1*id_pre)
v1_sin_pre1<-sin(2*pi*f_1*id_pre)
v2_sin_pre1<-sin(2*pi*f_2*id_pre)
pre_point<-data.frame(v1_cos=v1_cos_pre1,v1_sin=v1_sin_pre1,v2_sin=v2_sin_pre1)
hol.pred <- predict(tt,pre_point,interval = 'prediction', level = 0.95)
lcl_hol<-hol.pred[,'lwr']
ucl_hol<-hol.pred[,'upr']
#节假日范围=非节假日预测范围+节假日预测范围 ?
lcl_holtot<-lcl_ar[1]+lcl_hol
ucl_holtot<-ucl_ar[1]+ucl_hol
根据节假日模型预测节假日第一天数值的上下限,再加上非节假日模型预测的节假日第一天数值的上下限,两者总和就是节假日第一天预测值的上下限。模型也需要周期的检测,根据实际情况调参。