[deprecated] Rで大規模な実時間データのカテゴリー化
【2021-06-15 更新】本記事の内容は時代遅れです。 局所解なので、 アクセスログを見ながら優先度をつけてアップデートします。
実時間データの分析には色んなやっかいごとがあるように思えます。
例えば以下のデータは瞳孔経を300Hzで記録したもので、
時間と瞳孔経(mm)がそれぞれstamp列とpupil列にあります。
| subject | stamp | pupil | order | item | condition |
|---|---|---|---|---|---|
| P05 | 1 | 1.54 | 3 | 6 | d |
| P05 | 5 | 1.70 | 3 | 6 | d |
| P05 | 8 | 2.14 | 3 | 6 | d |
| P05 | 11 | 1.59 | 3 | 6 | d |
| P05 | 15 | 1.53 | 3 | 6 | d |
| P05 | 18 | 2.05 | 3 | 6 | d |
ただデータが多すぎるせいで、 可視化する時に重くなったり 偽陽性がでてしまったりします1。 あとはほとんど変化がないためデータとして冗長だったり、 機材によってサンプリング周波数が違ったり…といったところです。
もっと前向きなモチベとしては、 特定の区間にあるデータを同じ区間とみなし、 グループ全体で平均や和を計算したい、 なんてこともあるかもしれません。 しかしこのままだとどこが一つの区間かがわからないため不可能です。
そういう時に連続値を任意の間隔で区切ってカテゴリーを与える処理をビン分割なんて言ったりします。 そこで今回は使いまわせそうなビン分割関数を組んでみます2。
データの取得
まずはライブラリをimportします。 ただ以下のライブラリは書きやすくなる、読みやすくなるだけで実際は一つも要りません。
library(dplyr)
library(magrittr)
使うデータは大学のサーバーに置いてあるので
urlで取ってきてloadします3。
summaryで内容を確認できます。
最初に画像刺激を呈示してから10秒近いデータがとれてますね。
この被験者さんの瞳孔系は2mm前後というところみたいです。
# データ読み込み
load(url("http://phiz.c.u-tokyo.ac.jp/~kisiyama/gca2018/pupil.sample.Rdata"))
summary(pupil.sample)
# subject stamp pupil order item condition
# P05 :2909 Min. : 1 Min. :1.28 Min. :3 6:2909 a: 0
# P07 : 0 1st Qu.: 2584 1st Qu.:2.13 1st Qu.:3 b: 0
# P08 : 0 Median : 5154 Median :2.19 Median :3 c: 0
# P09 : 0 Mean : 5276 Mean :2.19 Mean :3 d:2909
# P10 : 0 3rd Qu.: 8017 3rd Qu.:2.27 3rd Qu.:3
# P11 : 0 Max. :10440 Max. :3.45 Max. :3
# (Other): 0
今回のミソはstamp とpupilです。
もう少し細かく見るためにheadして最初の20行をみます。
pupil.sample %>% head(20)
| subject | stamp | pupil | order | item | condition |
|---|---|---|---|---|---|
| P05 | 1 | 1.54 | 3 | 6 | d |
| P05 | 5 | 1.70 | 3 | 6 | d |
| P05 | 8 | 2.14 | 3 | 6 | d |
| P05 | 11 | 1.59 | 3 | 6 | d |
| P05 | 15 | 1.53 | 3 | 6 | d |
| P05 | 18 | 2.05 | 3 | 6 | d |
| P05 | 21 | 1.53 | 3 | 6 | d |
| P05 | 25 | 1.85 | 3 | 6 | d |
| P05 | 28 | 1.66 | 3 | 6 | d |
| P05 | 31 | 1.69 | 3 | 6 | d |
| P05 | 35 | 1.76 | 3 | 6 | d |
| P05 | 48 | 1.28 | 3 | 6 | d |
| P05 | 61 | 2.03 | 3 | 6 | d |
さて本題です。 今回の目標は stampを区切ってカテゴリを与えることです。
例えば20msごとに区切りたい場合は
以下のようにbin(三列目)を追加したいわけです。
20,40,60,80と20ms区切りのカテゴリがそれぞれの行に与えられています。
そして以下のようにビン分割するには hist関数が役立ちます。
| subject | stamp | bin | pupil | order | item | condition |
|---|---|---|---|---|---|---|
| P05 | 1 | 20 | 1.54 | 3 | 6 | d |
| P05 | 5 | 20 | 1.70 | 3 | 6 | d |
| P05 | 8 | 20 | 2.14 | 3 | 6 | d |
| P05 | 11 | 20 | 1.59 | 3 | 6 | d |
| P05 | 15 | 20 | 1.53 | 3 | 6 | d |
| P05 | 18 | 20 | 2.05 | 3 | 6 | d |
| P05 | 21 | 40 | 1.53 | 3 | 6 | d |
| P05 | 25 | 40 | 1.85 | 3 | 6 | d |
| P05 | 28 | 40 | 1.66 | 3 | 6 | d |
| P05 | 31 | 40 | 1.69 | 3 | 6 | d |
| P05 | 35 | 40 | 1.76 | 3 | 6 | d |
| P05 | 48 | 60 | 1.28 | 3 | 6 | d |
| P05 | 61 | 80 | 2.03 | 3 | 6 | d |
hist関数でビン分割
histはたぶんRを勉強し始める時にもっとも最初のころに触れる関数で
ヒストグラムをだしてくれます。
上のテーブルのstampが0msから60msまでの範囲を
20msでビン分割すると下のような図ができます。
ちなみに上の図は以下で描けます。
bin = 20
pupil.sample %>% filter(0<=stamp,stamp<=60) %$%
hist(stamp, breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin))
ここで大切なのは
一つのビンの中にある個数がそのまま列の数になるという点です。
図の0–20は6つ、21–40は5、41–60は1であることが分かりますが、
これが上の表のbinと対応しています。
このヒストグラムから上のビン分割した列を作るためには ステップが3つ必要です。
- 対応するビンの最大値(20,40,60,80…)を取得
- 0–20の個数、21–40の個数とそれぞれのビンにある要素の数を取得
- ビンの個数分、対応する最大値を持つ行を作成
histはヒストグラムを描くための情報を保持しているので、
それらの情報にアクセスして1と2を実現させます。
bin取得
まずは横軸の情報、つまり何を基準にhistがビン分割したか、
を取得できます。
ただそのままやると毎度ヒストグラムが生成されるので、
plot=FALSEと指定して図を出力しないようになります。
bin = 20
bins =
pupil.sample %>% filter(0<=stamp,stamp<=60) %$%
hist(stamp,
breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin),
plot=FALSE) %$%
breaks[-1]
bins
# [1] 20 40 60
結果は20, 40, 60と20msでビン分割した結果の数字が取れています。 あとはそれぞれのビンの頻度を求めます。
bin個数取得
各ビン内のデータが幾つあるか、という情報はcountsに入っています。
bin = 20
n_bins =
pupil.sample %>% filter(0<=stamp,stamp<=60) %$%
hist(stamp,
breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin),
plot=FALSE) %$%
counts
n_bins
# [1] 6 5 1
結果は6,5,1なので、0–20までが6つ、20–40までが5つ, そして40–60までが1つという結果が返ってきました。 あとはこれを反映したベクトルがあれば 先ほどの表を再現できます。
mapply
0–20までが6つ、20–40までが5つ, そして40–60までが1つ という結果を反映したベクトルを生成するためには ビンに対してビン個数を対応させたベクトルをつくります。 これに関しては一度見た方が分かりやすいと思います。
mapply(rep,bins,n_bins)
# [[1]]
# [1] 20 20 20 20 20 20
#
# [[2]]
# [1] 40 40 40 40 40
#
# [[3]]
# [1] 60
これはrepという関数をbinsとn_binsの各要素のそれぞれに
適応させています。つまり
rep(bins[1],n_bins[1]),rep(bins[2],n_bins[2]),...
というのを全ての要素分おこないます。
これをunlistするとベクトルになるので、
それを列に加えます。
最終的な関数と利用方法
以上をまとめた関数が以下になります。
# take a list of timestamp and bin, return the same number of rows with binned stamps
binner = function(stamp, bin){
bins = hist(stamp,
breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin),
plot=FALSE) %$%
breaks[-1]
n_bins= hist(stamp,
breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin),
plot=FALSE) %$%
counts
mapply(rep,bins,n_bins) %>% unlist %>%
return
}
後はこの関数をどうやって適用するかですね。
subject, itemにgroup_byしてあげると、
続く操作をsubjectとitemごとに行えます。
その操作というのがビン分割なので、stampとbinを引数に与えて
ビン分割したベクトルを返します。
そのままcbindするともともとのデータフレームと順序が違うので
(subject, itemをベースにした順序になっている)、
もとのデータフレームにarrangeを当ててあげています。
pupil.sample %>% group_by(subject,item) %>%
do(data.frame(bins=binner(.$stamp,bin=20))) %$%
cbind(bins,arrange(pupil.sample,subject,item)) %>% head(13)
する以下のような結果になります。
ちゃんとstampに対応したbinが作られていますね。
| bin | subject | stamp | pupil | order | item | condition |
|---|---|---|---|---|---|---|
| 20 | P05 | 1 | 1.54 | 3 | 6 | d |
| 20 | P05 | 5 | 1.70 | 3 | 6 | d |
| 20 | P05 | 8 | 2.14 | 3 | 6 | d |
| 20 | P05 | 11 | 1.59 | 3 | 6 | d |
| 20 | P05 | 15 | 1.53 | 3 | 6 | d |
| 20 | P05 | 18 | 2.05 | 3 | 6 | d |
| 40 | P05 | 21 | 1.53 | 3 | 6 | d |
| 40 | P05 | 25 | 1.85 | 3 | 6 | d |
| 40 | P05 | 28 | 1.66 | 3 | 6 | d |
| 40 | P05 | 31 | 1.69 | 3 | 6 | d |
| 40 | P05 | 35 | 1.76 | 3 | 6 | d |
| 60 | P05 | 48 | 1.28 | 3 | 6 | d |
| 80 | P05 | 61 | 2.03 | 3 | 6 | d |
まとめ
あとはsubject,item,condition,binあたりで
group_byしてsummariseすれば平均や標準偏差を求めるなど、
好きな操作ができます。
参考
- https://note.nkmk.me/python-pandas-cut-qcut-binning/
- https://www.ibm.com/support/knowledgecenter/ja/SSLVMB_24.0.0/spss/base/idh_bander_gating.html#idh_bander_gating