時間序列與數據挖掘
一、實驗說明
1. 環境登錄
無需密碼自動登錄，系統用戶名shiyanlou，密碼shiyanlou
2. 環境介紹
本實驗環境采用帶桌面的Ubuntu Linux環境，實驗中會用到：
1. LX終端（LXTerminal）: Linux命令行終端，打開后會進入Bash環境，可以使用Linux命令
2. GVim：非常好用的編輯器，最簡單的用法可以參考課程Vim編輯器
3. R:在命令行輸入‘R’進入交互式環境，下面的代碼都是在交互式環境運行
4. 數據：在命令行終端輸入以下命令：
# 下載數據 wget http://labfile.oss.aliyuncs.com/courses/360/synthetic_control.tar.gz # 解壓數據到當前文件夾 tar zxvf synthetic_control.tar.gz
3. 環境使用
使用R語言交互式環境輸入實驗所需的代碼及文件，使用LX終端（LXTerminal）運行所需命令進行操作。
完成實驗后可以點擊桌面上方的“實驗截圖”保存並分享實驗結果到微博，向好友展示自己的學習進度。實驗樓提供后台系統截圖，可以真實有效證明您已經完成了實驗。
實驗記錄頁面可以在“我的主頁”中查看，其中含有每次實驗的截圖及筆記，以及每次實驗的有效學習時間（指的是在實驗桌面內操作的時間，如果沒有操作，系統會記錄為發呆時間）。這些都是您學習的真實性證明。
二、課程介紹
1. R中的時間序列數據
2. 將時間序列分解為趨勢型、季節型以及隨機型
3. 在R中構建ARIMA模型，並用其預測未來數據
4. 介紹動態時間規整（DTW），然后使用DTW距離以及歐式距離處理時間序列從而實現層次聚類。
5. 關於時間序列分類的三個例子：一個是使用原始數據，一個是使用經過離散小波變換（DWT）后的數據，另外一個是KNN分類。
三、課程內容
1、R中的時間序列數據
使用ts這個類可以抽取等時間距離的樣本，當參數frequency=7的時候說明選取樣本的頻率單位是周，以此類推，頻率為12和4分別生成月度和季度數據。具體實現如下：
# 生成1-30的整數，頻率是12也就是月度數據，從2011年3月開始 > a <- ts(1:30, frequency=12, start=c(2011,3)) > print(a) # 轉換為字符串型 > str(a) # 輸出該時間序列的屬性 > attributes(a)
執行結果如下：
2、時間序列分解
時間序列分解就是將時間序列按趨勢性、季節性、周期性和不規則性依次分解。其中，趨勢部分代表的是長期趨勢，季節性指的是時間序列數據呈現季節性波動，周期性指的是指數據呈現周期性波動，不規則部分就是殘差。
下面講解一個關於時間序列分解的例子，使用的數據AirPassengers由Box & Jenkins國際航班從1949年到1960年的乘客數據組成。里面有144個觀測值。
> plot(AirPassengers) # 將數據預處理成月度數據 > apts <- ts(AirPassengers, frequency=12) # 使用函數decompose()分解時間序列 > f <- decompose(apts) # 季度數據 > f$figure > plot(f$figure, type="b", xaxt="n", xlab="") # 使用當前的時間區域給月份賦予名稱 > monthNames <- months(ISOdate(2011,1:12,1)) # 使用月份名稱標記x軸 # side=1表示設置x軸，at指的是范圍從10-12，las表示分割的單位刻度為2 > axis(1, at=1:12, labels=monthNames, las=2) > plot(f)
結果如下：
上圖中，'observed'標注的表表示原始的時間序列分布圖，往下第二個圖是顯示數據的呈現上升趨勢圖，第三個季節圖顯示數據受一定的季節因素影響，最后一個是在移除趨勢和季節影響因素后的圖。
思考：R里面還有哪些包，哪些函數可以實現時間序列的分解，並試着使用那些函數實現分解，並將分解結果進行對比。
3、時間序列預測
時間序列預測就是基於歷史數據預測未來事件。一個典型的例子就是基於股票的歷史數據預測它的開盤價。在時間序列預測中有兩個比較經典的模型分別是：自回歸移動平均模型（ARMA）和差分整合移動平均自回歸模型（ARIMA）。
下面使用單變量時間序列數據擬合ARIMA模型，並用該模型預測。
# 參數order由（p,d,q）組成，p=1指的是自回歸項數為1，list內容是季節seasonal參數
> fit <- arima(AirPassengers, order=c(1,0,0), list(order=c(2,1,0), period=12))
# 預測未來24個月的數據
> fore <- predict(fit, n.ahead=24)
# 95%的置信水平下的誤差范圍（U,L）
> U <- fore$pred + 2*fore$se
> L <- fore$pred - 2*fore$se
# col=c(1,2,4,4)表示線的顏色分別為黑色，紅色，藍色，藍色
# lty=c(1,1,2,2)中的1，2指連接點的先分別為實線和虛線
> ts.plot(AirPassengers, fore$pred, U, L, col=c(1,2,4,4), lty = c(1,1,2,2))
> legend("topleft", c("Actual", "Forecast", "Error Bounds (95% Confidence)"),col=c(1,2,4), lty=c(1,1,2))
預測結果如下：
4、時間序列聚類
時間序列聚類就是基於密度和距離將時間序列數據聚類，因此同一個類里面的時間序列具有相似性。有很多衡量距離和密度的指標比如歐式距離、曼哈頓距離、最大模、漢明距離、兩個向量之間的角度（內積）以及動態時間規划（DTW）距離。
4.1 動態時間規划距離
動態時間規划能夠找出兩個時間序列的最佳的對應關系，通過包'dtw'可以實現該算法。在包'dtw'中，函數dtw(x,y,...)計算動態時間規划並找出時間序列x與y之間最佳的對應關系。
代碼實現：
> library(dtw) # 平均生成100個在0-2*pi范圍的序列idx > idx <- seq(0, 2*pi, len=100) > a <- sin(idx) + runif(100)/10 > b <- cos(idx) > align <- dtw(a, b, step=asymmetricP1, keep=T) > dtwPlotTwoWay(align)
a與b這兩個序列的最佳對應關系如下圖所示：
4.2 合成控制圖的時間序列
合成控制圖時間序列數據集‘synthetic_control.data’存放於當前工作目錄'/home/shiyanlou'下，它包含600個合成控制圖數據，每一個合成控制圖都是由60個觀測值組成的時間序列，那些合成控制圖數據被分為6類：
- 1-100：正常型；
- 101-200：周期型；
- 201-300：上升趨；
- 301-400：下降趨勢；
- 401-500：上移；
- 401-600：下移。
首先，對數據進行預處理：
> sc <- read.table("synthetic_control.data", header=F, sep="")
# 顯示每一類數據的第一個樣本觀測值
> idx <- c(1,101,201,301,401,501)
> sample1 <- t(sc[idx,])
合成控制圖時間序列樣本數據分布如下：
4.3 使用歐式距離層次聚類
首先從上面的合成控制圖每一類數據中隨機選取10個樣本進行處理：
> set.seed(6218) > n <- 10 > s <- sample(1:100, n) > idx <- c(s, 100+s, 200+s, 300+s, 400+s, 500+s) > sample2 <- sc[idx,] > observedLabels <- rep(1:6, each=n) # 使用歐式距離層次聚類 > hc <- hclust(dist(sample2), method="average") > plot(hc, labels=observedLabels, main="") # 將聚類樹划分為6類 > rect.hclust(hc, k=6) > memb <- cutree(hc, k=6) > table(observedLabels, memb)
聚類結果與實際分類對比：
圖4.1
圖4.1中，第1類聚類正確，第2類聚類不是很好，有1個數據被划分到第1類，有2個數據划分到第3類，有1個數據划分到第4類，且上升趨勢（第3類）和上移（第5類）並不能很好的被區分，同理，下降趨勢（第4類）和下移（第6類）也沒有被很好的被識別。
4.4 使用DTW距離實現層次聚類
實現代碼如下：
> library(dtw) > distMatrix <- dist(sample2, method="DTW") > hc <- hclust(distMatrix, method="average") > plot(hc, labels=observedLabels, main="") > rect.hclust(hc, k=6) > memb <- cutree(hc, k=6) > table(observedLabels, memb)
聚類結果如下：
圖4.2
對比圖4.1和4.2可以發現，由后者的聚類效果比較好可看出在測量時間序列的相似性方面，DTW距離比歐式距離要好點。
5、時間序列分類
時間序列分類就是在一個已經標記好類別的時間序列的基礎上建立一個分類模型，然后使用這個模系去預測沒有被分類的時間序列。從時間序列中提取新的特征有助於提高分類模系的性能。提取特征的方法包括奇異值分解（SVD）、離散傅里葉變化（PFT）、離散小波變化（DWT）和分段聚集逼近（PAA）。
5.1 使用原始數據分類
我們使用包'party'中的函數ctree()來給原始時間序列數據分類。實現代碼如下：
# 給原始的數據集加入分類標簽classId > classId <- rep(as.character(1:6), each=100) > newSc <- data.frame(cbind(classId, sc)) > library(party) > ct <- ctree(classId ~ ., data=newSc,controls = ctree_control(minsplit=30, minbucket=10, maxdepth=5)) > pClassId <- predict(ct) > table(classId, pClassId) # 計算分類的准確率 > (sum(classId==pClassId)) / nrow(sc) > plot(ct, ip_args=list(pval=FALSE),ep_args=list(digits=0))
輸出決策樹：
5.2 提取特征分類
接下來，我們使用DWT（離散小波變化）的方法從時間序列中提取特征然后建立一個分類模型。小波變換能處理多樣的頻率數據，因此具有自適應性。
下面展示一個提取DWT系數的例子。離散小波變換在R中可以通過包'wavelets'實現。包里面的函數dwt()可以計算離散小波的系數，該函數中主要的3個參數X、filter和boundary分別是單變量或多變量的時間序列、使用的小波過濾方式以及分解的水平，函數返回的參數有W和V分別指離散小波系數以及尺度系數。原始時間序列可以通過函數idwt()逆離散小波重新獲得。
> library(party)
> library(wavelets)
> wtData <- NULL
# 遍歷所有時間序列
> for (i in 1:nrow(sc)) {
+ a <- t(sc[i,])
+ wt <- dwt(X=a, filter="haar", boundary="periodic")
+ wtData <- rbind(wtData, unlist(c(wt@W, wt@V[[wt@level]])))
+ }
> wtData <- as.data.frame(wtData)
> wtSc <- data.frame(cbind(classId, wtData))
# 使用DWT建立一個決策樹，control參數是對決策樹形狀大小的限制
> ct <- ctree(classId ~ ., data=wtSc,controls = ctree_control(minsplit=30, minbucket=10, maxdepth=5))
> pClassId <- predict(ct)
# 將真實分類與聚類后的類別進行對比
> table(classId, pClassId)
5.3 K-NN分類
K-NN算法也可以用於時間序列的分類。算法流程如下：
- 找出一個新對象的k個最近鄰居
- 然后觀察該區域內擁有相同屬性的數量最多的類代表該區域所屬的類
但是，這種直接尋找k個最近鄰居的算法時間復雜度為O(n**2)，其中n是數據的大小。因此，當處理較大數據的時候就需要一個有效索引的方式。包'RANN'提供一個快速的最近鄰搜索算法，能夠將算法的時間復雜度縮短為O(n log n)。
下面是在沒有索引的情況下實現K-NN算法：
> k <- 20 # 通過在第501個時間序列中添加噪聲來創建一個新的數據集 > newTS <- sc[501,] + runif(100)*15 # 使用‘DTW’方法計算新數據集與原始數據集之間的距離 > distances <- dist(newTS, sc, method="DTW") # 給距離升序排列 > s <- sort(as.vector(distances), index.return=TRUE) # s$ix[1:k]是排行在前20的距離，表哥輸出k個最近鄰居所屬的類 > table(classId[s$ix[1:k]])
輸出結果如下：
由上圖輸出表格數據可知，20個鄰居中有19個數據都屬於第6類，因此將該類時間序列划分到第六類時間序列中。
**思考**：經過學習這么多時間序列分類，請思考以上時間序列分類方法的利與弊。
更多數據挖掘教材請來實驗樓學習。