Fitbitから取得した心拍データで時系列の異常検知を試してみる
井出先生の「異常検知と変化検知」を読んで,自分でも試してみたいと思ったんですが,あいにくちょうどいい時系列データが手元にないなーと思ってました.そんな折,データサイエンスLT祭りの発表の中に,Fitbitデータを可視化するものがあって*1,これはちょうどいいということで試してみましたよというていのエントリになります.
- 作者: 井手剛,杉山将
- 出版社/メーカー: 講談社
- 発売日: 2015/08/08
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (2件) を見る
Fitbitってなによ
Fitbitが何かしらない人のために一応説明しておくと,最近はやりの活動量計です.私が持っているのは,心拍が取得できるタイプのやつです.風呂に入るとき以外は一日中つけっぱなしで,睡眠とか運動とかを自動で判定してくれるので,手間がかからず便利です.自転車に乗っている時間もちゃんと判定してくれるのは,判定アルゴリズムよくできているなーという感じですね*2.ちなみに運動や睡眠の判定は,データをスマホに転送するタイミングで一気に行われるので,このときに時系列の状態推定アルゴリズムを走らせてるんだろうなぁと思います.その一方で歩数に関しては,手を振ったりしてもほとんどカウントされず,ちゃんと歩いている場合にのみ歩数が増えるところをみると,バンドパスフィルタでも使って地面に足が着地する際のパルス応答的な変化を取ってるんじゃないかなーと推測します*3.
【日本正規代理店品】Fitbit ワイヤレス活動量計+心拍計リストバンド ChargeHR Large Black FB405BKL-JPN
- 出版社/メーカー: Fitbit
- 発売日: 2015/04/24
- メディア: エレクトロニクス
- この商品を含むブログ (7件) を見る
Fitbitのデータを取得する
もちろん,Fitbitの各種センサから得られる生データを取得することはできません.ですがFitbitは開発者向けのAPIを提供しており,ここからデータを取得できるようです.また,もっとベタにWebサイトをスクレイピングしちゃうやり方もあるみたいですね.この辺りをサポートしてくれるパッケージをざっと調べたところ,Rなら fitbitScraper,Pythonなら python-fitbit あたりでしょうか.あとはRなら {httr} パッケージとか使って,OAuthからのAPI直叩きという方法もありそうです.
今回は,普通に {fitbitScraper} を使って,心拍データを取得することにしました.得られるデータは5分間隔のため,1時間で12点,1日で288点になります.以下のようなコードで,適当な日付の心拍データ変動を取得してグラフ化してみました.(fitbitScraper) だけでなく,あとで2軸グラフを記述するために {devtools} 経由で {plotflow} パッケージも入れておきます.(fitbitScraper) は www.fitbit.com へのログインメアドとパスワードが必要なので,持っていない人は取得しておいてください.
install.packages("fitbitScraper") devtools::install_github("trinker/plotflow") library(fitbitScraper) library(plotflow) library(ggplot2) 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)
明け方の寝ている時間は心拍数が50弱/分で推移しており,起きてから犬の散歩にいったので100くらいまで上がっています.そのあとは通常の動きで,夕方前くらいにサイクリングにいったため,やはり100をコンスタントに超えています.そして24時に近づくにつれ心拍が落ち着いていくのが見て取れます.
こんな感じで,Fitbitのデータは割と簡単に取得できます.
特異スペクトル変換法による時系列データの異常検知
今回使うデータは,時系列の心拍データになります.ここで目的としたいのは,それまでの時間と比べて明らかにおかしな動きがみられたとき,つまり変化点を検知したいということです.井出先生の本にもいくつかの手法が解説されていますが,今回はその中で特異スペクトル変換法を取り上げたいと思います.
この手法は,ざっくりいうとある時点tを起点として,その前のw時点分のデータ(これを部分時系列データと呼びます)を1時点ずつ前にずらしていって,それをk個たばねたもの(仮にと呼びます)と,そこからL時点分後ろにずらして同じくk個たばねたもの(と呼びます)を比較します.両者が似ていれば,t時点で変化は起こっておらず,逆に両者が大きく異なっていたら,t時点付近で時系列データに変化が起こったと考えましょう,ということになります.
このは,1行につきw時点ぶんのデータが並び,それをk行分たばねたk x w行列になります.も同様です.知りたいのはとがどれだけ似ているか,ですが,その前段として, と のそれぞれについて,特異値分解で要約し,左特異値ベクトルをm次元分だけ取り出します*4.これにより,時系列の変化を考慮した上で,とのそれぞれを要約することができました.
そして得られた2つの特異値について,以下の計算をすることで,最終的な異常値を得ることができます.左辺の第2項は,両者が類似しているほど1に近い値になるので,これを1から引くことで非類似度,つまりこの場合はそれ以前からの変化の度合いを表す値が得られるわけです.
Rで特異スペクトル変換法を実装
さて,それでは実際に計算してみましょう.このあたりの式は,同じく井出先生の「入門機械学習による異常検知」を下敷きにしています.
- 作者: 井手剛
- 出版社/メーカー: コロナ社
- 発売日: 2015/02/19
- メディア: 単行本
- この商品を含むブログ (4件) を見る
# 定数 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="Anomary 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)
全部で288時点しないので*5,部分時系列の時点数は16くらいに絞っておきます.結果は以下の通りです.おおむね想定したとおりに,変化点が検知できているように思います,ただ後半のサイクリングの変化点が,微妙に後ろにずれちゃっているようにも思います.このあたり,細かいスパイクが悪い影響与えているかもしれないので,スムージングしてからスコアを出したほうがいいのかなぁとも思います.
パラメタのチューニング
上の例では,すべての定数を決めうちにしちゃっているので,いくつかパラメタを変えてどうスコアが変わるかを見ていきたいと思います.ここではwとmについて試してみます.
まずはwからみていきましょう.wを8, 16, 32の3パターンに変えて,スコアを出してみたのが下の図になります.青が8,赤が16,緑が32になります.みるとわかるように,wが増えるほどスコアが減少して,変化を検出しづらくなっています.部分時系列が増えるんだから,当然そうなりますよね.
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="Anomary 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)
同様に,mについても1, 2, 3と変えてみてみたのが,以下の図になります.青が1,赤が2,緑が3になります.こちらも同様に,次元数が少ないほど結果がピーキーになってますね.このあたり,元の時系列データの特性に合わせて,適切な次元数を決めてあげる必要がありそうです.
そんなわけで,ライフログデータでいろいろ試してみることができますね.面白い.コードはgistにあげてあるので,よろしければどうぞ.
*1:@millionsmileさんの「チャップが山粕を一撃する方法」という発表なんですが,特に資料は公開されてないみたいですね...
*2:とはいえ,感覚的には15分以上連続で運転するぐらいでないと,サイクリングと判定をしてくれないので,誤判定を防ぐためなんでしょうけど,割と保守的だなーとは思います.
*3:そもそも歩数カウントはオンラインでちゃんとカウントされていくため,時間&空間計算量をそこそこ必要とする時系列の状態推定は厳しそうに思います.なので単純なフィルタリングをしてるだけなんじゃないかなー
*4:ここで特異値分解で得られる左特異値ベクトルのL2ノルムをとると1になります
*5:この日のデータは欠損があるっぽくて,286時点しかありませんでしたが.
Ansible経由でDockerイメージを作ってみる
これまで仮想化とかクラウドとか,そんなにお仕事で触ってなかったこともあって割と放置気味だったのだけど,さすがに少しは使えないとねということでちょいと試してみましたというお話.
以前に依頼を受けてWebアプリを作ったことがあって,これを1年くらい前にAnsible化してました.で,これをDockerに乗せてElastic BeanstalkにDokcerコンテナとしてデプロイしようと思ったわけです.
DockerfileとAnsible playbook
最初はplaybookの内容を,全部Dokcerfileに書きかえようと思ったんですけど,これが思った以上に面倒くさい.ていうかDockerfileってAnsibleみたいな丁寧なセットアッププロセスや便利なラッパーがなくて,割とベタに RUN コマンドでスクリプトそのまま書かなきゃいけないの,かなり嫌です.このあたりって,他の人たちはどうやって解決しているんでしょう.細かい粒度でプロセス切り分けて,Dockerfileが肥大化しないようにしているから,そんなに問題にならないってことなんですかね.Ansibleのほうがラッパーなぶん手軽だし,依存関係とかも
ということで,前の資産を生かす形で,Ansible playbookを拡張することにしました.手順としては,まずDockerコンテナを立ちあげるタスクを追加します.そして立ち上げたコンテナに従来のplaybookを流して,あとはイメージを保存,そしてElastic Beanstalkにコンテナをデプロイ,です.
Ansible playbookのタスク追加
sites.yml を以下のように記述します.docker_host に対するタスクとして,docker タスクを追加しています.
--- - hosts: docker_host become: yes remote_user: docker tasks: - name: deploy centos container docker: image: centos:centos6 name: iroha ports: 80:80 expose: 80 tty: yes - hosts: container connection: docker roles: - nginx - php - iroha - mysql
このあたり,Dokcerとの連携については以下のエントリに詳しいです*1.
Docker machineのセットアップ
Dokcer for Macが出たので,これをインストールしておきましょう*2.いくつかセットアップしておかなければいけないことがありますが,このあたりは以下のエントリをみて適当にやってください.docker-machineでちゃんとホストVMをセットアップしておかないと,ちゃんとDockerコンテナを立てることはできません.
playbookの実行とイメージの保存
ここまで準備できたら,あとはplaybookを流せば,Webアプリの載ったDockerコンテナが立ち上がります.
ansible-playbook -i hosts site.yml
docker psで実際にコンテナができてるか確認して,imageを作成します.
$ docker ps CONTAINER ID IMAGE COMMAND CREATED STATUS PORTS NAMES 3985cf72cf58 centos:centos6 "/bin/bash" 13 days ago Up 13 days 0.0.0.0:80->80/tcp iroha b6ca8a876c82 jupyter/datascience-notebook "tini -- start-notebo" 13 days ago Up 13 days 0.0.0.0:8888->8888/tcp compassionate_joliot 3879d9ab028a tokyor/rstudio "/init" 2 weeks ago Up 2 weeks 1410/tcp, 0.0.0.0:8787->8787/tcp awesome_noether $ docker commit -m"initial commit for irohajirui-sho container" 3985cf72cf58 iroha:latest sha256:1455b4022f1977361d78d96a23aaca3a0c287d0f6a6ae14bc7c5607d04945e4a
このあとは
できあがったイメージをDocker hubにでもあげて,Elastic Beanstalkあたりでデプロイすれば,すらばしい継続的デリバリー体制の完成,めでたしめでたし.というところなんですが,現状まだElastic Beanstalkでちゃんと動作していない次第.うまくいったら別エントリにでもしますかね... ほんとうはMySQLをRDSに分離して,アプリケーションはDockerfileに切り出して,DockerfileでのElastic Beanstalkへのデプロイにしたほうが,より洗練されている感はあります.
Fundamental Models for Forecasting Elections から考える,日本の選挙における当選者予測方法
参議院選挙が終わったと思ったら今度は都知事選と,まさに選挙の夏ですね.選挙といえば,20時の投票〆切と同時に発表される,マスメディア各社の当選者予測が風物詩です.開票率0%で当確が続々打たれる様は,まさに統計学+社会調査の面目躍如という感じがしますね.各社の当選者予測については,詳細は必ずしも明らかにされていないものの,大枠については様々なところで述べられています.読みやすい記事としては,以下の記事なんかがよくまとまっているのではないでしょうか.
基本としては,RDDベースの社会調査をもとに,過去データに基づいて値の補正や当選確率予測モデルの当てはめを行うことで,高い精度の当確予測を行っている訳ですね.各社は定期的に政党支持率等の調査を行うとともに,選挙期間は頻度を上げて調査を行います.その中で,明確な差があると判断され,かつ出口調査でも傾向が変わらなければ,基本的には20時の時点で当確を打つことが可能になります.こうした選挙調査の手法については,日経リサーチの鈴木さんが書いた「選挙における調査と予測報道」や,朝日新聞社の松田さんが書いた「調査手法転換時の対応と判断 : 2000年総選挙と2001年参院選挙の事例(選挙とOR)」,さらにNHKの仁平さんの「選挙と出口調査」あたりにまとまっています.いずれも結構細かく述べられていて,興味深いです.
さて,こうした社会調査ベースの手法も興味深くはあるのですが,それ以外に当選者予測を行う方法はないのでしょうか.とりあえずCiNiiで日本の論文を探した限り,ほとんど存在しないというのが実情のようです*1 *2.唯一に近く引っかかったのは,「Twitterにおける候補者の情報拡散に着目した国政選挙当選者予測」くらいですが,流石にこれは高い精度を求めた手法という訳でもないので,参考にはならなさそうです.
ということで,もう少し広く探してみたところ,以下の論文が引っかかりました.ということで今回はこの論文の紹介をしたいと思います.
Fundamental Models for Forecasting Elections
論文概要
この論文,正確にはアメリカ経済学会のカンファレンスペーパーとして2013年に発表されているものになります*3.著者はGoogleとMicrosoft Research勤務の経済学者で,基本的には計量経済学の手法が使われています*4.手法としては,アメリカの大統領選,上院選,州知事選などに対して,3ヶ月以上前に手に入るようなマクロデータ等を用いてprobit回帰を行っています.
この論文内の先行研究レビューでは,主に選挙直近に手に入るデータを用いた計量分析や,予測市場に基づいた予測の手法等があげられています.これら手法と比べて,選挙の3ヶ月以上前に当選者を予測することができる,州ごとの予測を既存研究よりも高精度で予測することができる,というのがこの論文の主張になります.
手法
予測する変数は,候補者ごとの得票率と勝利確率の2つになります.これらを,過去の選挙データや経済・政治に関するマクロデータを用いて予測します.得票率については線形回帰,勝率予測はprobit回帰を用いています.各選挙について,州ごとに結果を予測することで,選挙回数に比してサンプル数を増やすことが可能になり,かつ各州の個別要因を考慮することが可能になります.
またこのモデルの特徴は,基本的に民主党と共和党のみを重視したものになっていることです.アメリカの場合は,2大政党制が長く続いており,この形でモデルにしてしまうのが妥当といえます.具体的には,大統領選を予測するモデルの中で,現職大統領が民主党か否かという独立変数を用いていたります(民主党なら1, 共和党なら-1となるダミー変数です).
また,民主党共和党以外の候補については,そのまま独立変数としてモデルに組み込むような荒技を使っています.例えば1968年の大統領選では,アメリカ独立党のジョージ・ウォレス候補がジョージア州・アラバマ州等で勝利を収めていますが,これを「Wallace's vote share in state in 1968 if Southern state and year is 1972」という独立変数として用いています.こうしたことが可能なのも,これらの事象があくまで例外的で,2大政党制が基本的には崩れていないからでしょう.
結果
ここでは,一番精度のよかった大統領選のモデルを示します.線形回帰による得票率予測モデルは,以下の通りです.変数には納得間のあるものがたくさん並んでおり,過去選挙の州ごとの投票率の偏差などは,日本でもそのまま当てはめることができそうにも思います.Wallace, Anderson, Perotあたりは,特別な選挙の影響排除の項目ですね.これらを含んだ結果の調整済決定係数は とそれなりに高い値を示しています.
また,結果について,横軸に民主党候補者の予測当選確率/得票率を,縦軸に実際の得票率を示したのが下の図です.たしかに,それなりに綺麗に予測できているといえそうですね.
同様に,州知事選や上院選についても,同様の結果がありますが,そちらは論文本体を参照してくださいということで省略します.
雑感
割と綺麗にまとめられた論文でした.さて当然ここで気になるのは,この結果が日本の選挙予測に応用可能なのか,という点です.私自身の感想としては,これだけだと結構厳しいなぁというものです.というのは,日本の選挙は政党の統合や新設が非常に多くあり,かつ政党間の細かい駆け引きや個別事情が入り組んでいます.今回の参院選でも,民主党と維新の会が合流して民進党になりましたが,この時点で既存の予測モデルをそのまま当てはめるのが困難になります.もちろん独立変数に頑張ってパラメタを入れ込むことも可能かもしれませんが,この形を取ると独立変数がそれだけで50個とかできちゃうように思います.
また時代とともに結構党の性格やあり方も変わっていくように思え,stableなモデルをどこまで当てはめてよいものか,結構悩ましいようにも感じます.かといって,動的線形モデルみたいなものを考えれば良いかというと,そもそも選挙という事象は数年に1回しか行われない,非常にサンプル数の限られるドメインです.そのため複雑なモデルを設定することは,それ自体が大きな不確実性を生む要因となってしまいます.そのためちまたで流行のディープラーニングなんかも,選挙に関しては使う機会がなさそうだなぁと思っていたりします.
ということで,詰まるところ既存のマスメディア各社の予測スタイルは,それはそれで最適解なのかもしれないなぁと感じる次第です.でもアカデミアの人たちには,もう少し予測モデル頑張ってほしいなぁと思ったり思わなかったり.
*1:「選挙 予測」みたいなキーワードで検索しても,ほとんど何も引っかからない状態です.
*2:なぜまず日本で探したかというと,選挙制度は各国で非常に大きく異なるため,海外の事例がそのまま日本に当てはまるとは考えにくいからです.
*3:2014年に Fundamental models for forecasting elections at the state level としてElectoral Studies誌にまとめられていますが,例のごとくElsevierの雑誌は有料なので,こちらは呼んでいません.
*4:全然関係ないですけど,アメリカのIT企業が経済学者を研究部門に抱えていたりするの,知を大事にしている感じがしてよいなーと個人的には思っています.
マーケットデザインと受入保留アルゴリズム
最近,ブログエントリを書くときの枕が読んだ本のことが多いですが,今回も御多分に洩れずであります*1.
現実の中でのマーケットデザイン
つい先日まで,以下の本を読んでました.マーケットデザインという分野を過分にしてしらなかったんですが,大元はゲーム理論の一種だったんですね.著者のロスさんらが,自身で証明した理論について現実問題に当てはめ,実際にうまく機能することを示した,その過程が描かれていて,非常に面白い本でした.臓器移植や公立学校の希望選択,新米医者の医局配属といった,いわゆるマッチングの問題が,周到に仕組みをデザインすることで,全体的な満足度が高い状態で解決された例をみると,なるほどなーという感想しかでないです.
Who Gets What (フー・ゲッツ・ホワット) ―マッチメイキングとマーケットデザインの新しい経済学
- 作者: アルビン・E・ロス,櫻井祐子
- 出版社/メーカー: 日本経済新聞出版社
- 発売日: 2016/03/19
- メディア: 単行本
- この商品を含むブログ (2件) を見る
この本の中では,主に公共に関わるマーケットのマッチング問題について,主に述べられていました.ですが私自身がこれを読もうと思ったきっかけは,新しいサービスやプラットフォームを構築するときには,マーケットデザインの考え方が割とうまく当てはめられるんじゃないかな,と思ったからです.典型的にはネットオークションは一種のマーケットですし,その流れで今ではメルカリのようなもう少しシンプルな個人売買のプラットフォームも同様だといえます.特にWebを介したB2CやB2B2Cのサービスは,割と多くがマッチング問題に関わるものではないかと思います*2.
そしてこの本の中でくどいくらいに繰り返し言われていたのは,アルゴリズム以外の部分,つまりはそのマーケットがどういう人たちで構成されていて,どういうセンシティブな問題があり,どういった政治的事情が絡んでいるか,といった細部まで踏み込んで緻密なデザインを行わなければ,マーケットはうまく機能しないということでした.さっきのヤフオクやメルカリなんかも,レビュー機能であったり決済代行機能であったり,それ以外にもさまざまな工夫や調整を行うことで,厚みのあるマーケット*3を構築できたのだろうなぁと思います.神は細部に宿るではないんですが,世の中そんなに簡単にスパッと割り切れるものではないということなんだなぁと,なんかしみじみ思ってしまいました.
よくあるマッチング問題
この本の中でメインで述べられていたのが,受入保留アルゴリズムというものです.本中の例だと,よりよい就職先を希望する医者の卵と,よりよい新人を採用したい医局側が,どうやればうまくマッチングするかという話が例として出されていました.なにもせず自由競争に任せてしまうと,他の医局に取られる前に有望そうな新人を先に取ってしまえということで,青田刈りのような状況になってしまいます.その青田刈りが行き着いた先には,まだ専門教育をちゃんと受ける前の段階の学生に対するアプローチ合戦が起こり,やりたいことや能力の評価がうまくできずマッチングミスが逆に増えてしまう,という結果が待っていました*4.
これを防止するために協定が結ばれ,一定期間より前には医局側が採用することができないというルールが決められたりもしましたが,時が経つにつれ内密に口利きがされたりという事例が相次ぎ,なし崩し的に元に戻ってしまう,という状態だったようです*5.このような現象が起こるのは,単純な希望順のマッチングを行うだけだと,一部の人気医局に希望が集中して,第1希望にあぶれてしまったが最後,第2希望も第3希望もすでに埋まってしまっており,希望順位の低いところにいかなければいけない制度のせいだったりします.そうすると,本当の希望順リストを書くよりも,確実にいけそうなところを第1希望にみたいなことが横行することになってしまうわけです.なので,この状態を解決するためには,本心の希望順リストを書いた上で,それが正しく機能するマッチングシステムを作らなければいけないのですね.
受入保留アルゴリズム
そこで出てくる受入保留アルゴリズムは,簡単にいうと,第1希望のマッチングを行った際に,マッチングしたペアの選択がその場で確定するわけではなく,仮のマッチングとして保留しておき,続く第2希望のマッチングでやってきた新たな希望者と合わせて,再度マッチングを行う,というプロセスをとります.つまりこの場合,最初のステップでたとえ第1希望でうまく滑り込めたとしても,より優秀な人が後のステップでその医局を希望してきたら,その人がマッチングすることになります.そして最初に滑り込めてた人は,第2希望の医局にマッチング希望を出すわけです.これを繰り返して,全員がマッチングしたら,その時に初めて結果が確定します.厳密なアルゴリズムについては,このあたりの講義資料でも読んでいただければと思います.
ということで,実際に実装して試してみました.上の例で示した医者-医局マッチングは,複数の対象を選択することができる応用的な問題なので,ここではもっとシンプルに1対1のマッチング問題を試しています.これは
安定結婚問題 - Wikipedia として知られている問題で,例としてお見合いパーティーにおける男性→女性の告白を考えます*6.実際に実行した結果も下に貼っていますが,3回ほど回して,6ステップ,12ステップ,22ステップでそれぞれ10対10のマッチングが安定状態に至りました*7.
そんな感じで,数理的な問題を実問題に落とし込めると楽しいだろうなぁと思った,そんな夜でした.この本はさすがに一般書すぎるので,もう少しメカニズムデザインについてちゃんと書いてある書籍を読みたいですね...何かおすすめあるんでしょうか.
*1:まぁ,お仕事に関する部分で書けるようなことがあまりない,というのが実際のところなのですが.
*2:Web以前からやられてますが,有名なリクルートのリボンモデルなんかは,まさにですよね.未読ですけど,須藤さんのこの記事なんかには,ちゃんとした分析が書いてあるのかなーと思います.
*3:この本の中では,多くのプレイヤーが集まり,活発に取引が行われている状況のことを,「厚みのある」という表現であらわしています
*4:過度な市場信奉者に対しての煽りでもないですが,シンプルな市場メカニズムに任せていれば最善の結果が得られるわけでは,必ずしもないということですね.
*5:なんか我が国の新卒就職市場を見ているようですね.
*6:フェミニズム的な意味で微妙な例だとは思いますが,他にパッとわかりやすい例が思い浮かばなかったので,ご勘弁ください.
*7:そもそも数学的に,このアルゴリズムは必ず安定状態に至ることが証明されています.ただ安定ではあるが,最善の結果とは限らない点はポイントです.正確には,voter側にとって最高,picker側にとって最悪の形で安定状態になります.
上司をマネジメントする (Managing your boss)
先日のETLの記事の中でも軽く触れたんですが,ホントの意味でデータが組織で活用されるためには,組織全体がデータを使って意思決定をする組織構成になっていないといけない,というのが最近強く感じることです.キンボール先生の本では,このあたりビジネスサイドのキーパーソンを必ず巻き込めと書かれているわけですが,逆にいうとそれしか書いてなかったりします.そんな折,TLでユウタロスさんが以下のようなツイートをしていました.
上司をマネジメントするという視点については、Managing Your Boss という論文で研究もされており、知人のハーバードビジネススクール出(本物)の経営者も在学中に読んでいたく感銘を受けたと言ってました。有効だと僕も思います。https://t.co/3gtnDiecJP
— ユウタロス Ver.1.3.4 (@bishop_ring) 2016年3月26日
データドリブンな組織を作るというのは,まさに組織のトップを抱き込んで動かしていくわけで,つまりは自分の上司やその上司,さらにその上司にまで首を縦に振らせていく必要があるわけです.部下のマネジメント法というのは,そこらの本屋にも結構置いてあるわけですけど,逆に自分の上司をマネジメントという観点でみるのは新鮮だなーと思って,GWの間に読んでみました*1.
Managing Your Boss (Harvard Business Review Classics)
- 作者: John J. Gabarro,John P. Kotter
- 出版社/メーカー: Harvard Business Review Press
- 発売日: 2008/01/08
- メディア: Kindle版
- この商品を含むブログを見る
ざっくりとしたレビュー
見出し的には以下のとおり.
- Misreading the boss-subordinate relationship
- Understanding the boss
- Understanding yourself
- Developing and managing the relationship
- Compatible work styles
- Mutual expectations
- A flow of information
- Dependability and honesty
- Good use of time and resources
この本を通しての一貫しているのは,上司と部下の関係は相互作用的な対人関係である,という社会心理学的な視点です*2.例えば成果を上げてのし上がってきた優秀なマネージャーが,気難しいボスとの関係がうまくいかず,結局事業が失敗してしまった場合に,気難しい上司にも問題はあるかもしれないが,このマネージャーは上司とうまく協調して物事を進めていくスキルが不足していた,という見方をします.とはいえ職場の力関係が示すように,基本的には部下が上司に対してどう付き合っていくかというのが,主要な観点になります.
Understanding the bossという見出しがあるように,まずは上司の性格,意思決定の仕方,ワークスタイルといった点を把握するところから始まります.その上で,上司がスムーズに業務をできるような形に自分のスタイルを調整して合わせていく,というのが基本的な流れです.そして面白いのは,これは功利的・政治的な問題というよりは*3,上司とうまく連携できることで,自分の仕事がうまく進められ,ひいては組織全体によい結果をもたらすという部分が強調されている点です.上司に媚びへつらうというよりも,自分の仕事をスムーズに進めるためには,そもそも上司と良好な関係を築けることが大事であるという話です*4.
もう少しだけ具体的な話でいうと,ワークスタイルを揃える(e.g. 文書で報告されるのが好きか,口頭が好きか),相互の期待を揃える(上司が何を期待しているかを把握する),情報の流し方を調整する(大事な問題だけを伝えるか,細かい情報も逐一伝えるか,ネガティブな情報を伝えるか)といった事柄が挙げられます.それ以外にも誠実であれ,とかどうでもいい話で上司の時間を浪費すると,本当に大事な案件に時間を割いてもらえない,とかそういった話が述べられています.
雑感
読んで思ったのは,これ自体は直属の上司との関係について書かれてはいるけど,組織内のキーパーソンとの関わり方全般について当てはまるよなぁ,ということです.組織にはいろいろな人がいるし,機能別組織であれば部署によって目指す目標は全く違うし,ある部署の目標と,3つ隣の部署の目標が正反対という場合だってあるわけです*5.なので複数の部署をまたいで物事を動かすには,相手の部署が何を目標としているか(何を期待しているか)を知った上で,自分のやりたいことが相手にとっても得になるように調整して,協調することが大事なわけです.これってまさにこの本で書かれていることですよね.
このことを敷衍していくと「情けは人のためならず」の心がけに行き着くかなーと思います.心理学的にはそれこそ古典のチャルディーニの「影響力の武器」なんかでも返報性の法則として書かれているように,こちらから何がしかの恩恵を施した場合,相手はどこかでそれをお返ししたくなるわけです.まぁ実際問題としては,そこまで打算的になってやるというより,一緒に仕事する相手が良い気持ちで仕事できる方が,自分も気分いいよねという単純な話でもいいかなと思ったりするわけですが.
- 作者: ロバート・B・チャルディーニ,社会行動研究会
- 出版社/メーカー: 誠信書房
- 発売日: 2014/07/10
- メディア: 単行本
- この商品を含むブログ (4件) を見る
あと,Harvard Business Review自体のまとめ記事がこちらにまとまっているので,合わせて読むとよいかなと思います.
ビッグデータの成熟期に改めて見直したいETL
Hadoopが出てきてから10年,ビッグデータという言葉が流行り始めてからでも5年以上が経ち,2016年現在では,Hadoopエコシステムを使ったデータ活用が当たり前のものとしてあります.とはいえ巷に出回っているビッグデータ活用事例というのは,綺麗な上澄みだけをすくい取っていたり,リリースしたてのピカピカのときに発表されていたり,というのが大半で,それが結構個人的に気に食わなかったりします.
ビッグデータが当たり前のものになっている現在においては,単に作っただけで価値があるというフェーズは過ぎ去っていて,継続的に運用しながら価値を生み出し続けることが,非常に重要な問題だと思います.特にビッグデータ界隈はミドルウェアやツールの陳腐化が激しく,またビジネス自体の変化速度も過去と比べてどんどん速くなっているわけで,そういった変化に対応していくためには,また別のスキルが必要とされるのではないでしょうか.このあたりについては,今年のHadoop conference JapanでされていたClouderaによるETLの概説,ドワンゴやDeNAの事例はまさにそういった問題を扱っています.みんな対外発表では綺麗なところを見せたいわけで,なかなかこうした実情が共有されることは少ないように感じています.
www.slideshare.net
www.slideshare.net
www.slideshare.net
でもこうした問題は,ビッグデータにより初めて起こったものではありません.データウェアハウス,ETLという言葉で,エンタープライズ業界では何十年も前から取り組まれており,ノウハウが積み上げられている分野でもあります*1.ビッグデータ活用自体がある程度成熟してきた今こそ,改めて長期的な運用を想定した,ETLの仕組み構築が注目されるべき時期に来ているのではないでしょうか*2.
データウェアハウスとETL
と,ちょっと挑発的な書き出しで始めましたが,要するに最近ラルフ・キンボールのThe Data Warehouse Toolkitを読んでいて,非常に感銘を受けたというお話です.このキンボール先生,スタースキーマで有名なディメンショナルモデル*3を提唱した人で,データウェアハウスの分野ではビル・インモンと並び最も有名な人の一人です.このThe Data Warehouse Toolkitは第3版まで出ていて,初版は日本語訳も出版されています*4.私は両方読んだんですが,第3版のほうがはるかに体系化されて,様々なノウハウが整理されているので,英語に苦がなければ原書で読むことをお勧めします.
The Data Warehouse Toolkit: The Definitive Guide to Dimensional Modeling
- 作者: Ralph Kimball,Margy Ross
- 出版社/メーカー: Wiley
- 発売日: 2013/07/01
- メディア: ペーパーバック
- この商品を含むブログを見る
- 作者: ラルフキンボール,Ralph Kimball,藤本康秀,岡田和美,下平学,伊藤磨瑳也,小畑喜一
- 出版社/メーカー: 日経BP社
- 発売日: 1998/05/14
- メディア: 単行本
- 購入: 2人 クリック: 68回
- この商品を含むブログ (8件) を見る
さて,では具体的な話に進みたいと思います.まずディメンショナルモデルとは,スタースキーマとはなんぞや,というお話ですが,このあたりの概説については id:EulerDijkstra さんの以下の記事が非常にわかりやすいので,一読をお勧めします.
ぶっちゃけビッグデータ時代になっても,基本的なデータ設計というのはディメンショナルモデルから大きく外れるものではありません.よくあるウェブサイトのアクセスログベースのユーザアクセス動向の分析についても,生のhttpdログをHDFSに突っ込んで,それを加工してユーザアクセスファクトテーブルを作ります.次にユーザマスタ,商品マスタを業務系のDBからSqoopで吸い出してきて,適当に加工整形してユーザディメンション,商品ディメンションテーブルを作ります.あとは,Hiveクエリでも書いてユーザセグメント毎の商品購入動向を集計するわけです.簡単ですね.
長期的な運用を前提とした仕組みづくり
ビッグデータ活用,みたいなプロジェクトが立ち上がる際,うちの会社でもログ解析ができるようにしたい,そこから何か有用な知見が得られるんじゃないか,みたいなフワッとした要件で始まることって割とあると思うんですけれども,これって結構つらい話です.キンボール先生が本の中で再三繰り返しているのは,データウェアハウス構築においては,ビジネス側の要件が一番大事で,それを達成するための仕組みづくりをしなければいけない,ということです.Hadoopがいくらコモディティサーバ使えばいいっていっても,またクラウドを使えば簡単に構築できるっていっても,それなりのデータ量でそれなりの分析をしようと思ったら,すぐに結構な金額が吹っ飛んでいきますし,費用対効果を考えないで立ち上げると,間違いなく長期的にドツボにはまっていきます.また,誰も使わないデータウェアハウスを作り続けることほど,エンジニアとして悲しいこともないですしね.
そんなわけで,プロジェクトをうまく進めるためには,必ずエンジニアサイドだけでなく,ビジネスサイドのキーパーソンを巻き込んで進めましょう,という話になります.その上で,ビジネスサイドの人たちが分析しやすいラベルをつけたディメンションテーブルを構築し,入り組んだ集計処理はあらかじめ行ったビューを作っておく,といった配慮が求められます.というのは,ややこしい集計作業をビジネスサイドの人が毎回正しく実行してくれると期待するのは,基本的に間違っていて,もし複数人の集計結果が異なっていたとしたら,問い合わせがくる先はデータウェアハウスの開発マネージャーのところだからです.そして正しい数字を担保しないと,せっかくデータウェアハウスを作っても,結局誰も使ってくれなくなってしまいます.
また,あらゆる基準は時とともに変化することを前提として,あらかじめ設計しておかなければいけません.特にディメンションテーブルについて,時間の経過とともにカラムの内容がゆるやかに変化していきます.これをキンボール先生はSlowly Changeng Dimension (SCD) と呼んで,本の中でも繰り返しその対処法を説明しています.例えばECサービスにおけるユーザ情報なんかがその典型で,ユーザの登録住所や電話番号なんて情報は,そう頻繁に変わるわけではないですが,数年とかのスパンではそれなりの変更があると考えられます.これをキチッとトレースできるようにしておかないと,過去に遡った分析を行うことができなくなります.
本の中ではこの問題に対して以下のような基本の4種類の対応法を示しています*5.
1. 上書きする
基本的に変更がなく,過去の情報を参照する必要がないような場合には,単純に上書きしてしまいます.
2. 新しい行を加える
これがベーシックかつ一番よく用いられる方法で,値が変わった場合には,別の行を追加します.ディメンションにはあらかじめ追加日と失効日の2カラムを持っておき,値が変わった場合に,既存の行に失効日を追記することで,その行が有効な期間を特定できるようにします.また最新フラグを示すカラムを追加しておくことで,簡単に最新状態の行だけを絞り込むことができるようにしておくと便利です.
3. 新しい属性を加える
最新の値と,その1個前の値を比べたい場合に,この手法を使います.やり方は簡単で,現在の値カラムと,直前の値カラムの2つをあらかじめ用意しておき,値が変更された場合には,既存の行を書き換えて対応します.
4. ミニディメンションを用いる
これは頻繁に変更が起こる場合に用いられます.ディメンションテーブルは,マスタとしてファクトテーブルと結合されるので,あまりにこのデータサイズが大きいと,パフォーマンスに影響が出ます.そのような場合には,変化する頻度の高いカラムだけを抜き出して,それら全ての組合せを網羅したディメンションテーブルをあらかじめ用意しておきます.それとは別に,変化する頻度の低いカラムは,従来のディメンションテーブルにまとめておきます.変更頻度によりディメンションテーブルをあらかじめ分割しておくことで,保守性が高い仕組みが構築できます.
またこの場合の結合キーですが,キンボール先生は,何の意味もない連番の数字か,ランダムハッシュを用いることを強く推奨しています.例えば商品ディメンションで考えるなら,もともと商品IDがあるのが普通なので,これを結合キーに使ってしまえば良いと考えがちですが,これだと商品IDに何がしかの変更が加わった際に,一気に窮地に陥ります.この手のナチュラルキーは,直感的である一方,サービスの統合やID体系の見直しといった,長期間の運用でたまに起こるイベントへの耐久性が皆無に近いわけです*6.
これら以外にも,ETLで必要とされるコンポーネントを34に分類して,その概要を述べていたり,典型的なシステム構成や事例に基づいたデータスキーマ,またデータウェアハウス構築プロジェクトの進め方といった多様な観点について述べられています.本のメイン部分は,Hadoop以前の時代に書かれているため,特にパフォーマンス的な部分では古いなーという部分があります.ただそれって,従来であれば3ヶ月ぶんしか扱えなかったデータが3年分扱えるようになった*7とか,集計軸を削減していたのがしなくてよくなったとかいう程度の違いしかなくて,データサイズが増える段階のどこかで,同じ問題にぶち当たります.いくらHadoopがスケールアウトするといっても,会社の予算は有限なので,サーバを無尽蔵に増やせるわけはありません.費用対効果のトレードオフを考えて,できる限りチューニングを行う必要があることには変わりがありません.
また,Hadoopで非構造化した生データを,加工整形する前に扱うことができるという点についても,使い所は絞らないといけないです.システム関連系では,生のデータをストリーム処理したり,ミニバッチ的に処理したりで,うまく活用してレコメンドだのに活かしたり,というのはあり得る使い方ですけど,それはあくまでデータサーエンティストやデータマイニングエンジニアみたいな,一部の専門家たちに限定されるべきです.特にビジネスサイドの人たちが触るデータについては,時間をかけて加工整形して,元データの変更部分をETLで吸収する必要があることに変わりはないし,そうしないと正しいデータ活用,データに基づいた意思決定にはつながっていかないんじゃないかなーと思っています.
そんなこんなで,とりとめなくETLについて書いてきましたが,The Data Warehouse Toolkitはとても良い本なので,みなさん読みましょう,ということが言いたかっただけのエントリでした.
*1:ただエンタープライズのこうしたノウハウって,あんまり世の中にシェアされることがないようには思っています.まぁB2Bの受託開発がメインだと,なかなかノウハウをシェアするメリットもありませんからね...
*2:ちなみにこのあたり,Treasure Dataとかのマネージドなデータウェアハウスサービスを使っていれば構築運用部分の悩みは解消されるにせよ,どうやってデータを持つ点は使う側が考えなければいけない以上,決してETLが必要ないわけではないです.
*3:コメントで,ディメンジョンでなくディメンションでは,といただいたので,修正しました.その通りですね,お恥ずかしい...
*4:すでに絶版ですが,Amazonなんかで中古を買うことはできます.表紙の古臭さが,時代の流れを感じさせますね.
*5:ほんとうは,それらを組み合わせた手法が3つ紹介されており,合計7つだったりはします
*6:このあたり,私自身も実務で実感しています.
*7:つまり,3年以上のデータを扱おうと思った場合には,結局パフォーマンス的な問題が現れてくる
optim()で階層構造を持った多変数の最適化をしてみる
最適化周りの処理について実装する必要が出たので,optim() を調べて使ってみましたよ,という話.
optim() って関数の形で表わせさえすれば,結構なんでも自由にできるっぽくて便利です.特に,階層的な構造を持ったデータの最適化もできるのは大きいです.
例
階層化ということで,ここで例に挙げるのは,比率の最適化です.イメージとしては,各クラスにおける学級委員長の選出を思い浮かべてもらえればと思います.
各クラスには複数の立候補者がいて,その支持率が既知のものとします*1.各クラスの生徒に対して,支持している候補者と,目指すクラスの方向性についてのアンケート聞いたとしてください.その得られたアンケート結果に対しなにがしかの係数をかけて,各立候補者の支持スコアを算出します.で,それを全立候補者のスコア合計で割ってあげれば,各立候補者の支持率が擬似的に出せます.
そうすると,クラスごとに算出した支持率と,実際の支持率の両者が得られます.アンケート結果にかける係数を調整することで,支持率の予測と実際の乖離を0に近づけましょう,ということになります.この場合,係数がかかるのは各クラスの生徒ですが,知りたいのはクラスごとの支持率という,クラス単位のパラメタだという階層構造になるわけです.
サンプルデータの作成
とりあえず,以下のようにクラスデータを作っておきます.わかりやすくするために,クラスは2つ,各クラスの候補者も2人としておきます.ここでポイントは,予め各候補者のカラムを作っておくことです.各生徒の支持対象はわかっているので,ここでは1または0の値で表しちゃいます*2.
library(dplyr) # 選択肢1, 2のデータを作成して結合 # 項目は以下の通り # クラス番号 # 候補者1支持ダミー # 候補者2支持ダミー # 価値観1 # 価値観2 location1 <- matrix(c(rep(1, 100), rep(1, 30), rep(0, 70), rep(0, 30), rep(1, 70), rep(1, 71), rep(0, 29), rep(0.3, 69), rep(0.7, 31)), 100, 5) location2 <- matrix(c(rep(2, 100), rep(1, 50), rep(0,50), rep(0, 50), rep(1,50), rep(0.1, 40), rep(0.8, 20), rep(0.5, 40), rep(0.2, 20), rep(1, 80)), 100, 5) d <- as.data.frame(rbind(location1, location2)) colnames(d) <- c("class_id", "person1", "person2", "value1", "value2") # クラス1, 2の支持率データ # 項目は以下の通り # クラス番号 # 候補者1支持率 # 候補者2支持率 s <- as.data.frame(t(matrix(c(1, 0.8, 0.2, 2, 0.6, 0.4), 3, 2))) colnames(s) <- c("class_id", "support1", "support2") # 結合してデータセットを作成 ds <- d %>% inner_join(s, by="class_id")
最適化対象の関数の作成
やりたいのは各候補者のスコアを出して,それを既知の支持率と比較することです.まずはわかりやすくするために,1クラス分のデータだけでやってみたいと思います.
# 1クラスぶんにデータを絞る ds1 <- ds %>% filter(class_id==1) # 最適化する関数 predict_rating <- function(b, d, print=FALSE) { d_tmp <- d d_tmp$predict1 <- b[1]*d_tmp$person1*b[3]*d_tmp$value1*b[4]*d_tmp$value2 d_tmp$predict2 <- b[2]*d_tmp$person2*b[3]*d_tmp$value1*b[4]*d_tmp$value2 # クラスidでグループ化して合計スコアを出した上で,比率に変換する d_tmp2 <- d_tmp %>% group_by(class_id, support1, support2) %>% summarize(., predict1=sum(predict1), predict2=sum(predict2)) total <- d_tmp2$predict1+d_tmp2$predict2 d_tmp2$predict1 <- d_tmp2$predict1/total d_tmp2$predict2 <- d_tmp2$predict2/total if (print) { print(d_tmp2) } # 最少化する関数は二乗誤差 sum((d_tmp2$support1-d_tmp2$predict1)^2+abs(d_tmp2$support2-d_tmp2$predict2)^2) }
optim() で実行
ここまで用意ができたら,実際にoptim() で実行します.この場合,最適化したい係数は非負で一定の範囲の実数に抑えたいとします.そのような指定ができる手法は,どうやら "L-BFGS-B*3" だけらしいので,これを選択します.各変数は0から10の範囲を取るという制約を与えます*4.各係数の初期値には,平均と分散が1の正規分布から乱数を生成し,負の値は正にするかたちで用いています.
res <- optim(abs(rnorm(4, 1, 1)), predict_rating, d=ds1, method="L-BFGS-B", lower=0, upper=10, control=list(maxit=10000)) res
これを実行すると,以下のように結果が得られ,実際に最適化がなされていることが見て取れます.
> res <- optim(abs(rnorm(4, 1, 1)), + predict_rating, + d=ds1, + method="L-BFGS-B", + lower=0, + upper=10, + control=list(maxit=10000)) > res $par [1] 3.2007261 0.5497433 0.5569501 0.7339547 $value [1] 1.603087e-14 $counts function gradient 10 10 $convergence [1] 0 $message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" > predict_rating(res$par, d=ds1, TRUE) Source: local data frame [1 x 5] Groups: class_id, support1 [?] class_id support1 support2 predict1 predict2 (dbl) (dbl) (dbl) (dbl) (dbl) 1 1 0.8 0.2 0.7999999 0.2000001 [1] 1.603087e-14
ここまで結果が得られたところで,最適化対象のクラス数を2にしてみましょう.ds1ではなくdsを使う,というだけですが.そうすると,見事に以下のように最適化された結果が得られています.もちろん2クラスに同じ係数を当てはめているため,クラス1については当てはまりが多少悪くなっています.それでもクラス1,2ともにそれなりによい当てはまりになっていることがわかるかと思います.
> res <- optim(abs(rnorm(4, 1, 1)), + predict_rating, + d=ds, + method="L-BFGS-B", + lower=0, + upper=10, + control=list(maxit=10000)) > res $par [1] 2.42237544 0.52544405 1.42063132 0.04326498 $value [1] 0.005155603 $counts function gradient 8 8 $convergence [1] 0 $message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" > predict_rating(res$par, d=ds, TRUE) Source: local data frame [2 x 5] Groups: class_id, support1 [?] class_id support1 support2 predict1 predict2 (dbl) (dbl) (dbl) (dbl) (dbl) 1 1 0.8 0.2 0.7600352 0.2399648 2 2 0.6 0.4 0.6313148 0.3686852 [1] 0.005155603