機器學習：異常檢測算法Seasonal Hybrid ESD及R語言實現

　　Twritters的異常檢測算法（Anomaly Detection）做的比較好，Seasonal Hybrid ESD算法是先用STL把序列分解，考察殘差項。假定這一項符合正態分布，然后就可以用Generalized ESD提取離群點。

　　目標是檢測出時間序列數據集的異常點，如圖所示，藍色線是時間序列數據集，紅色是圈是異常點。

　　

　　

 

 

 

 

 

 

 

 

 

 

　　R語言實現如下，一些依賴包需要install.packages("")或者手動在cran社區下載(注意依賴包的下載)。本人github下載源碼。

　　1 主函數是，包含了主要邏輯，加載數據集，IMSHESD算法檢測異常點和畫出數據集和異常點。IMSHESD算法是主要功能，下面詳細介紹。

library(zoo)
path_data="E:/Develop/Rstudio/IMS-H-ESDS/3151.csv"
path_sear="E:/Develop/Rstudio/IMS-H-ESDS/IS-H-ESD.R"

source(path_sear)
data <-read.table(path_data,sep=",",skip=1)
print(data)
test_data=IMSHESD(data)
plot(x=data[[1]], y=data[[2]],xlab="time",ylab="value",col="blue",type="l")
lines(x=test_data[[1]], y=test_data[[2]],col="red",type="p")

　　2 IMSHESD算法是主要邏輯如下，通過Fourier轉換自動求得時間序列的季節周期peri(必須滿足數據集的長度>2*peri+1才可以應用時間序列分析)，按照季節周期對數據集做划分，然后應用anmodetection異常檢測算法探測異常。

IMSHESD<-function(data,group_peri=60)
{
path_Fourier="E:/Develop/Rstudio/IMS-H-ESDS/Fourier.R"
path_data_group="E:/Develop/Rstudio/IMS-H-ESDS/data_group.R"
path_anmo_detection="E:/Develop/Rstudio/IMS-H-ESDS/anmo_detection.R"
source(path_Fourier)
source(path_data_group)
source(path_anmo_detection)
#data <-data_group(data,mode="median",group_period=60)
peri=Fourier_trans(data)
print(peri)
if(ncol(data)!=2)
  {
     print("The col of data must be two!")
     stop()
  }
if((2*peri+1)<length(data[[2]])){
   data_sep_length=ceiling(2*peri+1)
   print(data_sep_length)
   all_data <- vector(mode="list", length=ceiling(length(data[[1]])/(data_sep_length)))
   for(j in seq(1,length(data[[1]]), by=data_sep_length)){
      start_data <- data[[1]][j]
      end_data <- data[[1]][min(j + data_sep_length, length(data[[1]]))]
      if(j+data_sep_length<length(data[[1]])){
          all_data[[ceiling(j/(data_sep_length))]] <- subset(data, data[[1]] >= start_data & data[[1]] < end_data)
      }else{
         all_data[[ceiling(j/(data_sep_length))]] <- subset(data,data[[1]] >= data[[1]][length(data[[1]])-data_sep_length] & data[[1]] < end_data)
      }
   } 
   res=c()
   for(i in 1:length(all_data)) 
  {
      res_temp=anmodetection(all_data[[i]],anoms_per=0.1,period=peri,alpha=0.05)
      res=c(res,res_temp)
  }
  data_plot=rep(c(0),length(res))
  for(i in 1:length(res))
  {
     data_plot[i]=data[[2]][which(data[[1]]==res[i])]
  }
  anmo_point<-data.frame(res,data_plot)
}else{
    print("This is not a seasonal time series")
    stop()
}
}

　　3 Fourier轉換自動求得時間序列的季節周期peri。

Fourier_trans<-function(data)
{
   install.packages("TSA")
   library(TSA)
   p=periodogram(data[2])
   dd=data.frame(freq=p$freq,spec=p$spec)
   order=dd[order(-dd$spec),]
   top2=head(order,3)
   time=min(1.00/top2$f)
}

　　4 異常檢測主要邏輯anmodetection函數，需要規定異常點的上限10%，STL分解數據集：周期+趨勢+隨機噪聲=原始時間序列（分解方法有Twitters的Decompose和STL），殘差項根據正態分布（方差未知使用學生t分布），提取離散點。假設要檢測k個離群點，就對數據重復使用k次ESD檢驗，如果發現離群點就從數據里剔出去，然后在剩下的數據上重新檢測（Generalized ESD）。

anmodetection<-function(data,anoms_per=0.10,period=7,alpha=0.05,mode="addi")
{
   num <- length(data[[2]])
   #cat("num",num)
   num_anmo=trunc(anoms_per*length(data[[2]]))
   R_idx=rep(seq(0),length(data[[2]]))
   if(ncol(data)!=2)
     {
       print("The col of data must be two!")
       stop()
     }
   data_decompose <- stl(ts(data[[2L]], frequency =period),"periodic",robust = TRUE)
   seasonal_data <- data_decompose$time.series[,1]
   #cat("seasonal_data",seasonal_data)
   trend_data <- data_decompose$time.series[,2]
   random_decomp <- data_decompose$time.series[,3]
   data<-data.frame(times=data[[1]],count=(data[[2]]-seasonal_data-median(data[[2]])))
   func_med <- match.fun(median)
   func_mad <- match.fun(mad)
   numb_anom=0
   for(n in 1:num_anmo)
  {
     data_norm<-abs(data[[2]]-func_med(data[[2]]))
     data_mad<-func_mad(data[[2]])
     data_res<-data_norm/data_mad
     R_res<-max(data_res)
     max_temp_idx <- which(data_res == R_res)[1]
     R_idx_out=data[[1]][max_temp_idx]
     R_idx[n] <- R_idx_out
     data <- data[-which(data[[1]] == R_idx[n]), ]   
     p <- 1 - alpha/(2*(num-n+1)) 
     t <- qt(p,(num-n-1))
     thres <- t*(num-n) / sqrt((num-n-1+t**2)*(num-n+1))
     if(R_res>thres)
       {
        numb_anom <- n
       }
  }
  if(numb_anom>0)
  {
    R_idx <- R_idx[1:numb_anom]
  } 
  else 
  {
    R_idx = NULL
  }
  return(R_idx)
}

 

 　　實驗結果：紅色圈標出了異常點，異常點是10%的比例。

 

 

　參考

　　Twitters異常檢測方法，https://anomaly.io/blog/