mikutaifukuの雑記帳

個人的な雑記帳。データ分析とか読んだ本の感想とか。

Fitbitデータでスパース推定に入門してみる① 〜lassoで睡眠効率に影響を与える要因を探る〜

はじめに

2018年2月22日~2月23日(2日間)、統計数理研究所の2017年度リーディングDAT講座「【L-B2】 機械学習とデータサイエンスの現代的手法」(2017年度リーディングDAT)に参加してきました。(【L-B1】 実践ベイズモデリング」は抽選に外れたので来年あればリベンジしたい。) 従来の私ならこのような講座に参加するだけで満足しがちなのですが、しっかりと知識が定着するようにブログにまとめたいと思います。

講座内容は「カーネル法」「因果推論」「ニューラルネットワーク」「教師なし学習と行列分解」「スパース推定」と盛りだくさんでしたが、 今月末に阪大MMDS「機械学習・統計学 スプリングキャンプ」に参加するので、予習も兼ねてここでは「スパース推定(lassoを中心)」について取り上げます。

スパースとは

「スパース」とは「疎」という意味で、スカスカというイメージです。

私もこのような知的なツイートができるようになりたいです。

スパース推定とは

スパース推定とは統計学における推定方法の1つ(最小二乗法や最尤推定みたいな)で、特徴は

  • パラメータ推定
  • 変数選択(特徴量選択)

を同時に行うことにあります。データはものすごい大規模だけれど、意味ある情報はごく一部である時に使われます。具体例としては、以下のようなものがよく挙げられます。

  • 沢山ある遺伝子のうち、特定の病気に関係するものは少ない
  • 画像データでは、隣同士のピクセルの色は同じような色であることがほとんど

仕組みとしては、パラメータを0とすることで変数選択を可能とします。 f:id:mikutaifuku:20180304180916p:plain パラメータ推定値に0が含まれ、非0となる推定値がまばらにある様が、疎(スパース)のように見えるため、この推定方法をスパース推定とよびます。

一般的なパラメータ推定

f:id:mikutaifuku:20180302165636p:plain ただ、上記のような線形回帰モデルの場合、サンプルサイズより特徴量の方が大きい時(n < p)や、特徴量間の相関が高い時(多重共線性)は、最小二乗推定量は存在しないことが知られています。そこで、スパース推定の出番です。

スパース推定

スパース推定の考え方は非常にシンプルで上記に書いた最小二乗推定量を算出する式の後ろに、正則化項(罰則項)をつけて、パラメータが大きく推定されるのを防ぎます。特に、正則化項にL1ノルムを考えた場合をLassoと言います。Lassoはパラメータを0と推定することで、変数選択をすることが可能です。例えば相関が高い2つの説明変数があったとき、片方の説明変数のパラメータが0と推定されて、その変数はモデルから除外されます。「なぜLassoはパラメータを0に推定することができるのか?」はこの辺りの記事(https://research.preferred.jp/2014/05/group-lasso/)を参考にしていただければと思います。

f:id:mikutaifuku:20180302165849p:plain 医学領域では、DNAマイクロアレイデータ(遺伝子の発言量)やレセプトデータといった高次元小標本データは多く存在しており、スパース推定の応用が期待されています。

やってみた

タイトルにあるように、Fitbitデータに対してスパース推定を試してみました。fitbitデータ及び自作の生活行動データ(毎日Excelで手入力しています。)を利用して、睡眠効率に影響を与える要因を探ってみます。

データ

自作の生活行動データは欠損が非常に多く、データ数も多くありませんが気にせずやってみます。 項目は

  • totalTimeInBed:ベッドにいた時間(≠睡眠時間)
  • efficiency:睡眠効率(睡眠時間÷ベッドにいた時間)
  • Work :仕事の日は1
  • Drink:お酒を飲んだ日は1
  • Glass:眼鏡をかけていた日は1、コンタクトの日は0
  • Bath:湯船につかった日は1
  • Stretch:お風呂上がりにストレッチをした日は1
  • weekday:平日は1
  • step8000:8000歩以上歩いた日は1

となります。

> glimpse(df)
Observations: 43
Variables: 9
$ totalTimeInBed <int> 423, 357, 590, 426, 463, 426, 440, 417, 381, 588...
$ efficiency     <dbl> 0.9290780, 0.9607843, 0.9067797, 0.9436620, 0.95...
$ Work           <int> 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, ...
$ Drink          <int> 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, ...
$ Glass          <int> 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, ...
$ Bath           <int> 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, ...
$ Stretch        <int> 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, ...
$ weekday        <int> 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, ...
$ step8000       <int> 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, ...

本来であれば、しっかりと可視化をしてデータの分布などを確認すべきですが割愛します。

lasso

glmnetパッケージを使用すれば、簡単にlassoを実施する事ができます。

y <- unlist(select(df, efficiency))
x <- scale(as.matrix(select(df, -efficiency)))

library("glmnet")
fit <- glmnet(x, y, standardize = FALSE)
plot(fit, label = TRUE)

f:id:mikutaifuku:20180304132526p:plain

簡単にsolution path(解パス)を求める事ができました。solution pathとは、正則化パラメータλを変化させた時の、回帰係数の変化を示したものになります。パラメータが大きい時(罰則が大きい時)は、一部の変数のパラメータが0に推定されている事がわかります。

ただ、各線がどの変数を示しているのかが分からず、グラフもイマイチです。「ggplot2形式でsolution pathを描けないのか」と思い調べた結果、やっぱりありました。しかも綺麗にまとめてくださっていて素晴らしいですね。これがRの良いところだと思います。

qiita.com

library("ggfortify")
autoplot(fit)

f:id:mikutaifuku:20180304151133p:plain

裏側でggplotが動いているのだと思います。個人的には「背景白色派」なので、もしかして…と思い下記を実行しました。

autoplot(fit) + theme_pander() + scale_color_pander()

f:id:mikutaifuku:20180304151424p:plain

いい感じになりました。この辺りは雰囲気でggplot2をやっている人間なので、オブジェクトの考え方含め一度しっかりと勉強しようと思います。

ハイパーパラメータの選択

上記の結果を見てわかるように回帰係数は、ハイパーパラメータの選び方によって大きく変化します。従って、最適なハイパーパラメータを交差検証法(cross validation)等によって選択する必要があります。これも簡単にRで実行する事ができます。

set.seed(0)
cvfit <- cv.glmnet(x, y, standardize = FALSE, nfolds = 10)
autoplot(cvfit) + theme_pander() + scale_color_pander()

f:id:mikutaifuku:20180304152045p:plain

左側の点線が交差検証法で算出した2乗誤差を最小にするλ (= λmin) を表しています。
λminの時の回帰係数は以下のように取得にします。 weekday(平日)とstep8000(歩数が8000歩以上)が変数として選択されました。

> coef(cvfit, s = "lambda.min")
9 x 1 sparse Matrix of class "dgCMatrix"
                          1
(Intercept)    0.9431179145
totalTimeInBed .           
Work           .           
Drink          .           
Glass          .           
Bath           .           
Stretch        .           
weekday        0.0025008443
step8000       0.0006026092

結果の解釈

推定結果やsolution pathから自分なりに結果を解釈してみます。

  • weekdayは睡眠効率にプラス
    • 平日はベッドの中でダラダラする事が少ないからだと思われる
  • step8000は睡眠効率にプラス
    • よく運動している(歩いている日)はぐっすり眠りにつく事ができるので睡眠効率にプラスの影響を与えそう
  • Drink(飲酒)は睡眠効率にプラス
    • 家に帰ってきた瞬間に着替えてベッドに入り爆睡するから睡眠効率が良いのかも…。これは意外な結果でした。

λが小さい時は各変数間に相関関係が見られるので(多重共線性)、少し回帰係数の推定結果が微妙な気がします。Workとweekdayはほぼ同じ日フラグが立っているのですが、おそらく相関が強いのでプラスとマイナスに両極端な推定結果が見られているのだと思います(たぶん)。ただ、cross validationで得られた結果ではweekdayしか選択されなかったのでとりあえずは大丈夫なのかなと思います(私見)。

さいごに

今回のエントリでは、スパース推定入門という事で、Fitbitデータの睡眠効率に影響を与える要因をlassoで探索しました。この記事を作成するにあたって昨年参加した(C)スパース推定(平成29年度)の資料を参考に書きました。Rのコードの部分は岩波データサイエンス Vol.5の前半部分を参考にしました。今後この辺りの本をしっかり読んで(スパース推定法による統計モデリング (統計学One Point 6))理論面もしっかりと学ぼうと思います。

最後に個人的に気になった点、学びたい点をMEMOしておきます。

  • 説明変数が0か1のデータであった時にも、問題なくスパース推定出来るのか?説明変数の標準化について、ある行為を受けたか受けなったかみたいな01変数を正規化する事の意味・解釈の仕方がよくわからない。
  • ggplot2のオブジェクトについて諸々。雰囲気でRをやる人間なので、そろそろ真面目に勉強したい

分からないことだらけですね。まだまだINPUTが足らないようなので、引き続き勉強頑張ります。

なお、本エントリは雰囲気で書いた部分も多いので、間違いがあれば遠慮なくご指摘くださいませ。

mikutaifuku

Fitbitデータで因果を探索してみる 〜独立成分分析によるLiNGAMモデルの推定〜

はじめに

2018/2/11に日本経済新聞電子版にて以下のような記事が掲載されました。 f:id:mikutaifuku:20180218131146p:plain 元論文(日本経済研究センター JCER)の総論では、「博士増、生産性向上に結びつかず」と書かれており、大学教育や企業の活かし方に問題があるのでは、と問題提起をしています。

twitter等でも拡散され、盛り上がりました。 togetter.com

因果関係があるというのはなかなか難しく、疫学分野での因果推論の指標となっているヒルの因果関係ガイドラインを利用することが有用であると考えられています。
ヒルは事象Aが事象Bの原因であると結論付ける為に、以下の9つの規準

  1. 相関関係の強さ Aの生起とBの生起の間に強い相関関係がある
  2. 相関関係の一致性 相関関係の大きさは様々な状況で、対象や実証に利用よる手法が違っても一致している
  3. 相関関係の特異性 Bと「A以外に原因として想定される変数」の相関は高くない。またAと「B以外の結果変数」の相関も高くない
  4. 時間的な先行性 AはBに時間的に先行する
  5. 量・反応関係の成立 原因となる変数Aの値が大きくなると、単調に結果となる変数Bの値も大きくなる
  6. 妥当性 AがBの原因になっているという因果関係が生物学的に(または各分野の知見にもとづいて)もっともらしい
  7. 先行知見との整合性 これまでの先行研究や知見と首尾一貫している
  8. 実験による知見 動物実験等での実験研究による証拠がある
  9. 他の知見との類似性 すでに確立している別の因果関係と類似した関係・構造を有している

を満たしているかどうかチェックすることを推奨しています。*1
記事の内容が、上記のチェックリストにどれだけ満たしているかチェックしてみても面白いかもしれません。
また、今回の例では下記(ほんの一例)のように因果の向きは色々なパターンが考えられます。

f:id:mikutaifuku:20180218134133p:plain おそらく、上記論稿では過去の文献や調査などから、因果グラフが図でいうところの左上であるとみなして議論を進めたのだと思います。(たぶん)

このエントリでお話する「因果探索」とは、上記の中でどの因果グラフが正しいかをデータから推測することです。つまり、因果探索は因果グラフが未知の状態から、データを用いて仮説の候補を探索したりすることです。

ということで、例のごとく自分のFitbitデータを使ってやってみました。なお、本エントリは清水先生の書籍の1〜4章を参考にして書きました。

統計的因果探索 (機械学習プロフェッショナルシリーズ)

統計的因果探索 (機械学習プロフェッショナルシリーズ)

概要(やりたいこと)

Fitbitデータから取得できる、歩数・消費カロリーの因果関係を探ることです。
2017/01/01〜2017/12/31の1年間の歩数と消費カロリーデータをプロットしたのが以下の図です(横軸:歩数、縦軸:消費カロリー、相関係数:0.8)。 f:id:mikutaifuku:20180218030110p:plain
相関関係があるのは明らかですが、因果関係が以下のパターンのどれになるかを推測するのが目的です。

f:id:mikutaifuku:20180218135621p:plain
なお、今回は話を簡単にするために、

  • 歩数と消費カロリーに影響を与える他の未観測共通原因は存在しない

ことにします。未観測共通原因がある場合の話も今後したいと思います。

普通に考えれば「歩数を増やすと消費カロリーが増える」、つまり一番左の絵の「歩数→消費カロリー」が成り立ちそうですが、LiNGAMモデルを使って、これが本当に推測できるかをどうかを確かめてみます。

LiNGAMモデルとは

LiNGAMモデル(linear non-Gaussian acyclic model, LiNGAM)は因果グラフが識別可能なモデルです。ここで、識別可能とは

  • 因果グラフの構造が異なれば、観測変数の分布が必ず異なる。つまり、因果グラフが一意に推定できる。

ことを言います。LiNGAMは識別可能なモデルで、具体的には以下の条件を満たすモデルのことを言います。*2 f:id:mikutaifuku:20180218160655p:plain
今回の例は2変数(歩数・消費カロリー)なので、以下のLiNGAMモデルになります。係数行列のBを推定できれば、因果グラフを描くことができます。

f:id:mikutaifuku:20180218160714p:plain

独立成分分析によるアプローチ

計算の流れ

  1. x = Bx + e → (I - B)x = e → e = Wx より、独立成分分析を適用して、データ行列Xから復元行列Wを推定します

  2. 1で推定された復元行列 \hat{W}_{ICA}は、行の順序と尺度を一意に推測することができないので、適切な置換行列Pと尺度を変換する対角行列Dを推定します

  3.  {\hat{W} = \hat{D}^{-1} \hat{P} \hat{W}_{ICA}}で復元行列Wを推定します

  4.  {\hat{B} = I - \hat{W}}で行列Bを推定します

  5. 行列Bが厳密に下三角になるように(有向非巡回になるように)因果的順序を求めます

  6. 因果的順序に従い、係数行列Bを推定します

以上が、独立成分分析に基づくLiNGAMモデルの推定の流れになりますが、ここで書いたのはかなり雑な記述なので、詳細は書籍を読んでいただければと思います。

使用データ

上の散布図で示した通り、2017/01/01〜2017/12/31の1年間の歩数と消費カロリーデータを使用しました。欠損はありません。

図:歩数ヒストグラム f:id:mikutaifuku:20180218164353p:plain
図:消費カロリーヒストグラム f:id:mikutaifuku:20180218164450p:plain

歩数と消費カロリーのヒストグラムをみてわかるように、どちらもガウス分布に従ってそうで、LiNGAMモデルの仮定に当てはまらなさそうですが、書籍には

ただ、因果的順序を推測してしまえば、変数の推定には非ガウス性を利用する必要はありません。そのため、誤差変数の分布をガウス分布で近似することにして情報量規準を利用しても、実際の性能はそれほど悪くないでしょう

と書いているので、ここでは気にせずに進めます。*3

結果

pcalgパッケージ

独立成分分析でLiNGAMモデルを推定するにあたって、pcalgパッケージを使用しました。
graph, BiocGenerics, RBGLといったパッケージも同時にインストールする必要がありますが、こちらはCRANには登録されていないので、以下のサイトから圧縮ファイルをダウンロードして、Rにインストールしました。

www.bioconductor.org

データ

> tail(X)
       steps caloriesOut
[360,]  2747        2431
[361,] 12634        3290
[362,] 11148        3170
[363,]  3534        2424
[364,]  6282        2575
[365,]  6952        2734

推定結果

> library("pcalg")
> lingam(scale(X))
$Bpruned
          [,1] [,2]
[1,] 0.0000000    0
[2,] 0.8593356    0

となりました。$Bprunedは行列Bのことで、今回の例では  b_{12} =0となったので、「歩数→消費カロリー」の因果グラフを推定することができました。*4

さいごに

今回のエントリでは、Fitbitの歩数・消費カロリーデータに対して、因果グラフをLiNGAM法によって推定しました。
本来であれば、他の共通原因(心拍数etc)を加味する必要があるとは思いますが、次回以降に持ち越したいと思います。

最後に個人的に気になった点をMEMOしておきます。

  • 今回の分析は「1個人の反復測定データ」だったが、そもそもこのように因果探索してもいいのか
  • 連続変数でない時(0/1の二値変数etc)はできないのか
  • 結局、非ガウスでも因果的順序を使えば、LiNGAMモデルを推定できるのか
  • パラメータを推定できても、それが因果関係あるかは別じゃないのか

分からないことだらけですね。まだまだINPUTが足らないようなので、引き続き勉強頑張ります。

なお、本エントリは雰囲気で書いた部分も多いので、間違いがあれば遠慮なくご指摘くださいませ。

mikutaifuku

*1:調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)を引用しました。

*2:統計的因果探索 (機械学習プロフェッショナルシリーズ)P.95を参考に書きました。

*3:個人的なブログなので許してください。私の理解が浅いのもあり「因果的順序の決定→行列Bを推定」の場合は、非ガウス性の仮定がいるのかどうかわかりませんでした。

*4:行列Xをscaleで標準化しなかった場合は、「消費カロリー→歩数」となりました。係数を0に推定する過程で、標準化しておかないと不都合があるのでしょうか。この辺りも私にはよく分かりませんでした。

Fitbitの心拍数データで異常検知 〜自己回帰モデルによるアプローチ〜

概要

前回の記事に引き続き「入門 機械学習による異常検知―Rによる実践ガイド」を参考に、自分のFitbitの心拍数データを使って7章の時系列データの異常検知をやってみました。前回は特異スペクトル変換法を用いましたが、今回は自己回帰モデルを使用してやってみます。似たような記事は他にも色々あるので少し参考にさせていただきました*1。簡単ではありますが自分なりにまとめてみます。

使用データについて

今までの記事と同様に、Fitbitの睡眠時心拍数データを使用しています。

詳細は以下を参照してください。

mikutaifuku.hatenablog.com

自己回帰モデルとは

「現在の観測値が、過去の複数の観測値(例えばr個)の一次結合に従っている」と考えることを自己回帰モデルと言います。パラメータrのことを次数と言います。次数はAIC赤池情報量基準)という、モデル選択基準を使用します。詳細は以下で述べます。

計算イメージ
  • 下図のように、各時刻の観測値からデータ集合Dを作り、線形回帰と同様のアプローチで未知のパラメータを算出します。詳細な導出過程は書籍や上で紹介した記事を参考にしていただければと思います。
  • 長さTの時系列データからは、t = r + 1からTまでの、合計T - r個の観測値に対する予測式を作ることができます。
  • 多変量の場合もパラメータが行列、観測値がベクトルになりますが、基本的な考え方は同じです。

f:id:mikutaifuku:20180208221620p:plain

次数rの決定
  • モデルの複雑さは次数rで決まります。クロスバリデーションなどの手法でrを決めることができないので(時系列データだから)、モデル選択基準としてAIC赤池情報量基準)を用います。考え方や詳細は他の書籍*2に譲りますが以下のようになります。

    AIC(r) = -2 × (モデルの最大対数尤度) + 2 × (モデルのパラメータ数)

  • このAICが最小になるrを選びます。
異常度の定義
  • 算出したモデルから得られた予測値が、実測値と乖離していた時に異常があったという考え方に基づいて、異常度を以下のように定義します。

   f:id:mikutaifuku:20180208224359p:plain

  • 下図のように、予測値に対して実測値が乖離している場合は、異常度が高くなります。

f:id:mikutaifuku:20180208225200p:plain

結果

変化度の計算結果を以下に示します。綺麗に異常部位が検出されました。

f:id:mikutaifuku:20180208225559p:plain

ただ、実測値と予測値を比較してわかるように、あまりにもフィットしているので(過学習なのかな?)、少し心配です…。どこかミスっているのかな?

 

f:id:mikutaifuku:20180208225829p:plain

AICにより選択された次数r = 7だったので、過去7時点のデータを使用して、現在の値を予測することができます。

このような単純なデータであれば簡単な計算だけで、リアルタイム性も担保しながら異常度求めることができるので、個人的には結構いいように思います(私見)。

さいごに

今回のエントリでは、Fitbitの心拍数データに対して、自己回帰モデルによるアプローチで変化点を検出しました。この手法はわかりやすくていいですね。次は状態空間モデルによる異常検知をやってみようと思います。

ただ、少しワンパターンな気がしてきたので、次回は個人的に興味がある「因果探索」や「スパース推定」のお話をブログに書こうと思います。

 

 

※ 何か間違いがございましたらご指摘ください。

mikutaifuku

Fitbitの心拍数データで異常検知 〜特異スペクトル変換法による変化点検知〜

概要

前回の記事に引き続き「入門 機械学習による異常検知―Rによる実践ガイド」を参考に、自分のFitbitの心拍数データを使って7章の時系列データの異常検知をやってみました。前回は近傍法を用いましたが、今回は特異スペクトル変換法で変化点検知をやってみます。似たような記事は他にも色々あり*1、二番煎じではありますが自分なりにまとめてみます。

使用データについて

前回の記事と同様に、Fitbitの睡眠時心拍数データを使用しています。

詳細は以下を参照してください。

mikutaifuku.hatenablog.com

変化点検知とは

前回は外れ値検出でしたが、今回は変化点検知という問題になります。前回は以下の図でいうと赤◯の部分を検出することが目的でしたが、今回は青◯の部分を検出することが目的となります(多分この認識であっているはず…)。

f:id:mikutaifuku:20180207175241p:plain

特異スペクトル変換法による変化点検知

計算イメージ
  • 各時刻の観測値をそれぞれ扱うのではなく、w個の隣接した観測値をw次元のベクトルとして表します。
  • 時刻tの周りに、過去側と現在側において、k本の部分時系列データを使って、履歴行列X1とテスト行列X2を作ります。
  • 行列X1、X2に対して特異値分解を行い、上位m個の左特異ベクトルの行列を求めます。
  • 得られた行列から時刻tにおける変化度を計算します。

ここでは計算の流れを書いただけですので、詳細は書籍や上で紹介した記事を参考にしていただければと思います。ここの内容を理解するには基本的な線形代数の知識が必要です。

イメージ図(参考:「入門 機械学習による異常検知―Rによる実践ガイド」P.200)

f:id:mikutaifuku:20180207181243p:plain

具体化したイメージ図

f:id:mikutaifuku:20180207181307p:plain

図を見てわかる通り、時刻tの変化度を求める為には、t + L - 1までのデータが必要となります。したがって、よりリアルタイム性を求めるのであれば、窓幅wを小さくしたり、ラグLを小さくしたりする等、各パラメータを調整する必要があります。

結果

変化度の計算結果を以下に示します。前回の近傍法と比較して、とても綺麗に異常部位が検出されました。ただ、よく見ると変化点を若干遅れて検知していることがわかります。これはパラメータを調整することで改善できる部分もありますが、アルゴリズムの都合上限界がありそうです。  

f:id:mikutaifuku:20180207205043p:plain

計算できる範囲やリアルタイム性を考えると、一番下のw = 5, m = 2, k = 3, L = 2のパターンが良さそうです(私見)。

書籍にも書いてあるように特異スペクトル変換法は

  • 計算コストが高い
  • 確率分布からの意味づけがはっきりしない

という課題が挙げられています。ということで、これらを踏まえて次回は、同じデータに対して「自己回帰モデルによる異常検知」を実践してみようと思います。

さいごに

今回のエントリでは、Fitbitの心拍数データに対して、特異スペクトル変換法によるアプローチで変化点を検出しました。

RやPythonのプログラムも公開したいところですが、人様に見せるのが恥ずかしいレベルなので今回も諦めました。そのうち、コードを綺麗に書き直してアップしたいと思います。

今までは本を読むだけで満足していましたが、ブログにまとめることで思考が整理させれ、頭に定着しやすいような気がします。今後もサボらずに定期的に記事を書いていきたいです。

 

※ 何か間違いがございましたらご指摘ください。

mikutaifuku

Fitbitの心拍数データで異常検知 〜近傍法で異常部位検出〜

概要

「入門 機械学習による異常検知―Rによる実践ガイド」を読んで、自分でもやってみたいと思い、自分のFitbitの心拍数データを使って7章の時系列データの異常検知をやってみました。似たような記事は他にも色々あり*1、二番煎じではありますが自分なりにまとめてみます。

 Fitbitとは

www.fitbit.com

歩数・睡眠状態・心拍数等の様々な活動量を記録してくれるウェアラブルバイスです。一部の企業では、従業員にFitbitを配布して、健康管理なんかもしているみたいです。ちなみに私はFitbit Alta HRを使用していますが、最近発売されたfitbit ionicと呼ばれるスマートウォッチも気になっています。

ダッシュボートは以下のような感じです。

もちろんスマートフォンでも見ることができます。

f:id:mikutaifuku:20180205213157p:plain

使用データ

fitbitは公式APIを提供しているので、そちらからデータを取得しました。

以下の記事を参考にしました。(結構苦労しました)

qiita.com

睡眠時以外は、Fitbitを外していてデータが欠損していたり、運動して急激に心拍数が上がったりと少し扱いにくかったので、睡眠時の心拍数データに絞りました。

可視化

とりあえず以下のように1/8~1/14の睡眠時心拍数データを可視化してみました。

目的は、以下の図の赤丸で囲んでいるような異常値を、近傍法で検出することです。

f:id:mikutaifuku:20180205214759p:plain

  • fitbitの睡眠データから、睡眠開始/終了時間を取得できるので、その時間帯の心拍数のみ取得しています
  • 以下のように「睡眠ステージ」データも取得できますが、今回は使用しません。睡眠開始時間からの数分は「目が覚めた状態」のため心拍数が高くなっていると思われます。

    f:id:mikutaifuku:20180205220146j:plain

  • Rのggplot2を使用しています(赤丸は除く)

それにしても、社会人とは思えないほどよく寝ていますね。大学生のような睡眠時間です。(この週は休暇で休んだり、出張だったりで少しイレギュラーでした。)

近傍法による異常部位検出

計算イメージ
  • 各時刻の観測値をそれぞれ扱うのではなく、w個の隣接した観測値をw次元のベクトルとして表します。
  • 訓練用の時系列Dtrの各要素と、検証用の時系列Dのx(t)の距離を計算します。

イメージ図(参考:「入門 機械学習による異常検知―Rによる実践ガイド」P.196)

f:id:mikutaifuku:20180206152440p:plain

  • 上記で求めた距離(今回はユークリッド距離)のうち、最小なものk個選び、score(異常度)を算出します。scoreはk = 1の時は最近傍までの距離、k > 1の時は近傍距離の平均を考えることとします。

結果

"2018-01-11"を訓練用の時系列データとし、"2018-01-12"を検証用の時系列データとしました。

窓幅 w = 5, 10の2パターン、近傍数 k = 1, 10の2パターンで試した結果、以下のようになりました。

f:id:mikutaifuku:20180206144111p:plain

一番上の心拍数データで示した異常部位に対して対応した箇所でscoreが大きくなっていることがわかります。w = 5, k = 10でやった結果が一番良さそうです。

書籍にも書いてあるように、近傍法を用いた方法はノイズに弱いらしいので、もう少し別の方法を検討する必要があるそうです。 ということで、次回は同じデータに対して「特異スペクトル変換法」で、異常検知をしてみようと思います。

さいごに

今回のエントリでは、Fitbitの心拍数データに対して、近傍法によるアプローチで異常部位を検出しました。

みなさんが、ブログで書いているようにRやPythonのプログラムも公開したいところですが、人様に見せるのが恥ずかしいレベルなので今回は諦めしました。

インターネットから取得できるよく分からないデータよりも、自分のデータを使って分析をした方が楽しいですね。

 

mikutaifuku