Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Anomaly detection for Fitbit
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment