Anomaly detection for Fitbit

fitbit_anomaly.R

install.packages("fitbitScraper")
devtools::install_github("trinker/plotflow")
library(fitbitScraper)
library(ggplot2)
library(plotflow)
library(dplyr)

# Fitbitからのデータ取得
cookie = login(email="www.fitbit.com ログイン用メールアドレス", password="www.fitbit.com ログイン用パスワード")
hr_data <- get_intraday_data(cookie, what="heart-rate", date="2016-10-06")
names(hr_data) <- c("time", "heart_rate") # 変数名にハイフンが入っているのでリネーム
hr_data <- hr_data %>% filter(heart_rate > 0) # Fitbitを取り外しているときは，値が0になるため取り除く

hr_plot <- ggplot(data=hr_data, aes(x=time, y=heart_rate)) +
  geom_line(size=0.8) +
  ylim(min=0, max=120) +
  labs(x="Time", y="Heart rate") +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=16, face="bold"))
plot(hr_plot)


# 特異スペクトル変換法による異常値判定
w <- 16 # 部分時系列の要素数
m <- 2 # 類似度算出に使用するSVDの次元数
k <- w / 2 # SVD算出に用いる部分時系列数
L <- k / 2 # 類似度を比較する2つの部分時系列群間のラグ

get_anomaly_score <- function(d, w, m, k, L) {
  anomaly_score <- rep(0, length(d))
  for(t in (w+k):(length(d)-L+1)) {
    start <- t - w - k + 1
    end <- t - 1
    X1 <- t(embed(d[start:end], w))[w:1,] # t以前の部分時系列群
    X2 <- t(embed(d[(start+L):(end+L)], w))[w:1,] # 異常度算出対象の部分時系列群
    
    U1 <- svd(X1)$u[,1:m] # X1にSVDを行い上位m次元を抜き出す
    U2 <- svd(X2)$u[,1:m] # X2にSVDを行い上位m次元を抜き出す
    sig1 <- svd(t(U1) %*% U2)$d[1] # U1とU2の最大特異値を取得
    anomaly_score[t] <- 1 - sig1^2 # 最大特異値の2ノルムを1から引くことで異常値を得る
  }
  anomaly_score
}

hr <- hr_data$heart_rate
anomaly_score <- get_anomaly_score(hr, w, m, k, L)
as_data <- data.frame(time=hr_data$time, anomaly_score=anomaly_score)
as_plot <- ggplot(as_data, aes(x=time, y=anomaly_score)) +
  geom_line(color="red", size=0.8) +
  labs(title="", x="Time index", y="anomaly score") +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=16, face="bold"))
combined_plot <- ggdual_axis(lhs=hr_plot, rhs=as_plot)
plot(combined_plot)


# 特異スペクトル変換のパラメタによる変動
as_data_w08 <- data.frame(time=hr_data$time, anomaly_score=get_anomaly_score(hr, 8, m, k, L))
as_data_w16 <- data.frame(time=hr_data$time, anomaly_score=get_anomaly_score(hr, 16, m, k, L))
as_data_w32 <- data.frame(time=hr_data$time, anomaly_score=get_anomaly_score(hr, 32, m, k, L))
as_plot_w <- ggplot(as_data_w08, aes(x=time, y=anomaly_score)) +
  geom_line(color="blue", size=0.8) +
  geom_line(data=as_data_w16, aes(x=time, y=anomaly_score), color="red", size=0.8) +
  geom_line(data=as_data_w32, aes(x=time, y=anomaly_score), color="green", size=0.8) +
  labs(title="", x="Time index", y="anomaly score") +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=16, face="bold"))
combined_plot_w <- ggdual_axis(lhs=hr_plot, rhs=as_plot_w)
plot(combined_plot_w)

as_data_m1 <- data.frame(time=hr_data$time, anomaly_score=get_anomaly_score(hr, w, 1, k, L))
as_data_m2 <- data.frame(time=hr_data$time, anomaly_score=get_anomaly_score(hr, w, 2, k, L))
as_data_m3 <- data.frame(time=hr_data$time, anomaly_score=get_anomaly_score(hr, w, 3, k, L))
as_plot_m <- ggplot(as_data_m1, aes(x=time, y=anomaly_score)) +
  geom_line(color="blue", size=0.8) +
  geom_line(data=as_data_m2, aes(x=time, y=anomaly_score), color="red", size=0.8) +
  geom_line(data=as_data_m3, aes(x=time, y=anomaly_score), color="green", size=0.8) +
  labs(title="", x="Time index", y="anomaly score") +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=16, face="bold"))
combined_plot_m <- ggdual_axis(lhs=hr_plot, rhs=as_plot_m)
plot(combined_plot_m)