结合NCS的k-means对数据进行分类

https://towardsdatascience.com/using-b-splines-and-k-means-to-cluster-time-series-16468f588ea6

这篇文章中介绍了一种简化分类时间序列数据的方法，流程就是对数据使用NCS（自然三次样条插值），根据插值参数得到每组数据的系数，这个系数的大小可以自行决定，本文中使用的为10，但由于将数据做了去均值处理，所以最后得到参数大小为9。这里的操作非常类似对时间序列进行特征降维，对降维后的特征再进行k-means分类。

这样的好处就在于如果直接对时间序列数据进行分类，由于时间序列非常大，特别容易造成过拟合。并且，由于一些非监督学习模型的特点，分类模型可能只专注于时间序列上的值的大小，而忽略了时间序列间的上下移动与趋势，这就导致了分类模型存在误差。

使用NCS对数据进行“降维”，在使用k-means进行分类，这就是这篇文章的主要思路

First, load the data provided by UCI. It’s a 811-by-52 dataframe

library(splines)
library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.4     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(dbplyr)

## 
## 载入程辑包：'dbplyr'

## The following objects are masked from 'package:dplyr':
## 
##     ident, sql

#set.seed(4358) 
mydata = read.csv("https://archive.ics.uci.edu/ml/machine-learning-databases/00396/Sales_Transactions_Dataset_Weekly.csv")
ts = as.matrix(mydata[, 2:53])
row.names(ts) = mydata[,1]
ts = ts - rowMeans(ts)

Then I calculate the corresponding coefficients of each rows, and save it as matrix B.

B can be obtained by the formula below

\[B^t=(F^tF)^{-1}F^tX^t\] where F is 52-by-9 design matrix without the intercept, it’s a base matrix of NCS.

#F=as.matrix(rep(0,52*9),nrow=52,ncol=9)

x = seq(0, 1, length.out = ncol(ts))
F=ns(1:52,df = 9, intercept = FALSE)
F=t(t(F) - colMeans(F))
B=solve(t(F)%*%F)%*%t(F)%*%t(ts)
B=t(B)

mykm1=kmeans(B,6)

Plot the six clusters and its centerline

myK = 6
par(mfrow=c(2,3))
mycenters1=matrix(rep(0,myK*ncol(ts)),ncol=ncol(ts))
for(k in 1:myK){
  id=which(mykm1$cluster==k)
  plot(NA, xlim = c(1, ncol(ts)), ylim = range(ts), xlab = "Weeks", ylab = "Weekly Sales")
  B_i= B[id,]
  temp_mycenters1= colMeans(B_i)
  mycenters1[k,]=F %*% t(t(temp_mycenters1))
  for(i in 1:length(id)){
    lines(1:ncol(ts), ts[id[i],] , col="gray")
  }
  lines(1:ncol(ts), mycenters1[k,], col="red")
}

For comparison, we plot the 6 clusters generated by using original data

mykm2=kmeans(ts,6)
myK = 6
par(mfrow=c(2,3))
mycenters2=matrix(rep(0,myK*ncol(ts)),ncol=ncol(ts))
for(k in 1:myK){
  id=which(mykm2$cluster==k)
  plot(NA, xlim = c(1, ncol(ts)), ylim = range(ts), xlab = "Weeks", ylab = "Weekly Sales")
  ts_i=ts[id,]
  mycenters2[k,]= colMeans(ts_i)
  for(i in 1:length(id)){
    lines(1:ncol(ts), ts[id[i],] , col="gray")
  }
  lines(1:ncol(ts), mycenters2[k,], col="red")
}