Fitbitデータでスパース推定に入門してみる③ 〜心拍数データに対してfused lassoを使い、変化点を検出する〜
はじめに
前回の記事では、個々の変数を選択するための様々なスパース推定の手法をまとめ、Fitbitデータに適用することでそれぞれの手法を比較しました。
今回のエントリでは、少し趣向を変えて心拍数データに対して、fused lassoと呼ばれる隣合う変数の関係性を加味した手法使ってみます。なお、今回のエントリは、書籍の第３章のfused lassoの部分を自分なりに整理してまとめただけなので、丁寧に学びたい方は、書籍を読んでいただければと思います。
スパース推定法による統計モデリング (統計学One Point)
- 作者: 川野秀一,松井秀俊,廣瀬慧
- 出版社/メーカー: 共立出版
- 発売日: 2018/03/08
- メディア: 単行本
- この商品を含むブログを見る
fused lassoの概要
fused lassoは以下のような性質をもつデータに対して提案された手法です。
- p個の説明変数に順序関係がある
- 隣接する説明変数のいくつかは目的変数に同程度の寄与がある
最小二乗推定量を算出する式の後ろに、正則化項にL1ノルムと隣接する係数の差の絶対値を考えたものです。
この式を最小にすることでパラメータβを推定します。正則化項として係数の差の絶対値を考えているため、推定値が等しくなるように働くという仕組みです。
時系列データへの応用
Fitbitの睡眠時心拍数データ*1を使って、時間的に変化する値の推定を考えます。隣り合う時点の値はほぼ等しく、ときどきある時点でジャンプしています。
時刻tでの観測値をとし、
とモデリングします。εは観測ノイズです。隣り合う時点でβの値が変わりにくいという性質を、正則化項
で表現します。この項の効果で不連続に飛ぶ点の数を少なくしながら推定を行います。つまり、以下のような式を最小にする事で推定パラメータを求めます。
これは、上記の(１)式においてn = p, Xが単位行列、λ1=0を考えたのと同じです。*2
Rでやってみた
genlassoパッケージのfusedlasso1dを用いることでfused lassoによる推定を行えます。
データ
> glimpse(df) Observations: 414 Variables: 2 $ time <dttm> 2018-01-11 00:27:00, 2018-01-11 00:28:00, 2018-01-11 0... $ value <int> 102, 110, 106, 84, 79, 80, 79, 78, 77, 80, 77, 83, 83, ...
fused lasso
library(genlasso) y <- unlist(df$value) res <- fusedlasso1d(y) # fused lasso plot(res, col="grey") # plot
正則化パラメータλを変えた時のフィッティングを折れ線で示しています。L1正則化により変化点が検出できそうです。が、デフォルトのplotがわかりにくいので自作します。
resには様々なλ（全部で383個）において推定されたパラメータβが格納されています。適当なλを選択し（今回は50, 75, 100, 125, 150番目を使いました）、可視化しやすいようにデータを縦持ち（tidy化）にします。
library("dplyr") library("tidyr") step <- c(50,75,100,125,150) beta_df <- data.frame(res$beta)[, step] # 指定した列のβを抽出 beta_df$time <- df$time # 時点の列を追加 beta_df$HeartRate_ob <- df$value # 実測値の列を追加 colnames(beta_df) <- c(res$lambda[step], "time", "HeartRate_ob") # ヘッダー名をλの値に変更 # tidy化 plot_df <- beta_df %>% gather(key="lambda", value = "HeartRate_es", -time, -HeartRate_ob) %>% mutate(lambda = paste0("λ = ", round(as.double(lambda), 2)))
データが良い感じになりました。
> glimpse(plot_df) Observations: 2,070 Variables: 4 $ time <dttm> 2018-01-11 00:27:00, 2018-01-11 00:28:00, 2018-... $ HeartRate_ob <int> 102, 110, 106, 84, 79, 80, 79, 78, 77, 80, 77, 8... $ lambda <chr> "λ = 13", "λ = 13", "λ = 13", "λ = 13", "λ = 13"... $ HeartRate_es <dbl> 101.66667, 101.66667, 101.66667, 84.00000, 82.11...
可視化
library("ggplot2") library("ggthemes") ggplot(plot_df, aes(time, HeartRate_ob)) + geom_point(color = "darkgrey", alpha = 0.5) + geom_line(aes(time, HeartRate_es, color=lambda)) + facet_wrap(~lambda, scales = "free", ncol = 1) + theme_pander() + scale_color_pander() + theme(legend.position = 'none') + xlab("") + ylab("Heart Rate")
そこそこ綺麗なグラフが完成しました。λが大きいと直線に近づき、λが小さいと観測値に沿った折れ線が得られていることがわかります。λ=4のグラフを見ると4:00や6:00頃の変化点をしっかりと検出できています。
さいごに
今回のエントリでは、Fitbitの心拍数データに対してfused lassoを適用し、変化点の検出を試みました。ただ、実務においては、このように全てのデータが揃ってから変化点を検出しても遅く、要望として例えば「異常を５分前に予兆したい」みたいなことが挙げられます。が、この手法では残念ながら解決できません（おそらく）。という事で、これら課題をどうやったら解決できそうか、うまく問題設定を考えてチャレンジしてみようと思います。
このブログではパッケージを駆使して「やってみた」系の内容がほとんどです（てか全部）。自分でしっかり頭を使ってモデリングをしている訳でもなければ、自分で数式をしっかり理解した上でアルゴリズムを実装している訳でもありません。この状況は個人的にはなかなかヤバイと思いながらも、時間コストを考えるとこれが限界な気がしなくもありません。色々難しいですね、ちょっと色々考えてみます。
引き続き頑張ります。
mikutaifuku
