モモノキ＆ナノネと一緒にR言語でchangepointパッケージを使った変化点検出をやってみよう
changepointパッケージを使った簡単なデータの変化点検出
ナノネ、今回はR言語で変化点検出をやってみよう。
データが変化する箇所を検出するの？
うん。changepointパッケージを使うとデータの変化点を簡単に検出できるから試してみよう。
モモノキ、データは何に使う？いつものirisは向いて無さそうかも？
そうだね．．．まずは、乱数を使った疑似データで変化点検出ができるか試してみよう。
set.seed(1)
# テストデータ１
data1 <- c(rnorm(100, mean=0, sd=1),rnorm(100, mean=5, sd=1))
plot(data1, type='l')
title('data1')
前半100個は平均0で標準偏差1（分散1）。後半100個が平均5で標準偏差1（分散1）のデータだね。
テストデータが準備できたので、実際にchangepointパッケージを使ってみよう。まずはパッケージのインポート（changepointは必要に応じてインストールしてね）。
# install.packages('changepoint')
library(changepoint)
changekpointパッケージで変化点検出を実行にするには、主に3つの関数あるよ。
・cpt.mean() 平均だけの変化を検出
・cpt.var() 分散だけの変化を検出
・cpt.menavar() 平均と分散の変化を検出
用途が関数の名前通りだから、覚えられそう。
作成したテストデータ１は、平均が途中で変わっているだけだから、cpt.mean()が使えそうだね。もし、変化する点が1個だけと分かっている場合は、関数の引数にデータを渡すだけで変化点検出ができるよ。
cpt_mean_data1 <- cpt.mean(data1)
cpt_mean_data1
なんかデータが出力されたけど．．．うまく検出できたの？
関数の実行で返ってきたオブジェクトをplot関数に渡すと、変化点付きのグラフを自動で描いてくれるよ。
plot(cpt_mean_data1)
title('cpt_mean(data1)')
うまく検出できてる！しかも簡単。
検出しやすい疑似データだからね。実際のデータだとオプションパラメータを少しいじるケースも多いかも。設定できるパラメータが沢山あるので、詳しくはドキュメントなどを調べてみてね。
# テストデータ１をcpt.var()で試しに実行
cpt_var_data1 <- cpt.var(data1)
plot(cpt_var_data1)
title('cpt_var(data1)')
同じデータで、cpt.var()を実行してみた。期待通りに変化点検出されずOK。
# テストデータ１をcpt.meanvar()で試しに実行
cpt_meanvar_data1 <- cpt.meanvar(data1)
plot(cpt_meanvar_data1)
title('cpt_meanvar(data1)')
cpt.meanvar()なら、平均の変化もうまく検出できてる。
分散が途中で変化するテストデータ２を作ったから、同じように確認してみて。
set.seed(1)
# テストデータ２
data2 <- c(rnorm(100, mean=0, sd=1),rnorm(100, mean=0, sd=5))
plot(data2, type='l')
title('data2')
前半100個は平均0で標準偏差1（分散1）。後半100個は平均0で標準偏差5（分散25）のデータだね。
par(mfrow=c(3, 1))
# cpt.mean()
cpt_mean_data2 <- cpt.mean(data2)
plot(cpt_mean_data2, cpt.col='blue')
title('cpt.mean(data2)')
# cpt.var()
cpt_var_data2 <- cpt.var(data2)
plot(cpt_var_data2, cpt.col='red', cpt.width=3)
title('cpt.var(data2)')
# cpt.meanvar()
cpt_meanvar_data2 <- cpt.meanvar(data2)
plot(cpt_meanvar_data2, cpt.col='green')
title('cpt.meanvar(data2)')
par(mfrow=c(1, 1))
テストデータ２は分散だけ変化させているから、cpt.var()を使うと変化点をうまく検出できるみたい。
検出された変化点の場所を数値で知りたいときは？
変化点が何番目のデータか知りたいときは、cpts()関数に検出結果オブジェクトを渡せばOK。また、同じようにparam.est()関数を使うと変化点で分けたときの平均や分散も確認できるよ。
# 変化点の確認（データn番目）
cpts(cpt_var_data2)
# 平均、分散の確認
param.est(cpt_var_data2)
変化点が複数ある場合はどうするの？必ずしも1点だけとは限らないよね。
変化点が複数ある場合も使う関数は同じでいいんだけど、オプション引数を少し変更する必要があるよ。
まず最初に、変化点が複数あるデータを適当に準備しておこう。
set.seed(1)
# テストデータ３
data3 <- c(
rnorm(100, mean=0, sd=1),
rnorm(100, mean=5, sd=2),
rnorm(100, mean=9, sd=1),
rnorm(100, mean=3, sd=1)
)
plot(data3, type='l')
title('data3')
複数の変化点を検出したい場合は、関数の引数methodに"PELT"を指定すると複数検出にも対応できるよ。（methodは計算方法の指定で、PELT以外にもいくつか選択できるんだけど、詳しい内容はドキュメントを調べてね）
cpt_meanvar_data3 <- cpt.meanvar(data3, method='PELT')
plot(cpt_meanvar_data3, cpt.width=3)
title("cpt.meanvar(data3, method='PELET')")
複数の変化点も検出できてる！分散だけ変更したデータも作って、変化点検出を試してみる。
set.seed(1)
# テストデータ４
data4 <- c(
rnorm(100, mean=0, sd=1),
rnorm(100, mean=0, sd=2),
rnorm(100, mean=0, sd=7),
rnorm(100, mean=0, sd=9)
)
par(mfrow=c(2, 1))
plot(data4, type='l')
title('data4')
cpt_var_data4 <- cpt.var(data4, method='PELT')
plot(cpt_var_data4, cpt.width=3)
title("cpt.var(data4, method='PELET')")
par(mfrow=c(1, 1))
あれ、変化点の検出箇所が1個足りない？分散の差が小さすぎると検出は無理みたい。見た目でも違いが分かりずらいし。
検出力にも限界があるよ。データによっては、パラメータをうまく設定することで対処できるケースもあるので、機会があったらいろいろ試してみてね。
次はデータセット使って、changepointを試してみよう。今回使うのは時系列データで有名なNile。
たしかナイル川の年別水位だっけ？
Nile
Nileデータで変化点検出をやってみる。
cpt_nile <- cpt.meanvar(Nile, method='PELT')
plot(cpt_nile)
title('Nile')
# 変化点の確認（データn番目）
print(cpts(cpt_nile)) # インデックス
print((1871:1970)[cpts(cpt_nile)]) # 年
# 平均、分散の確認
param.est(cpt_nile)
meanvar(method='PELT')関数を実行したら、変化点を3つ自動検出したみたい。見た目だと1個でもよさそうなんだけどなぁ。
変化点が1つだけと決まっているなら、methodにamoc（または指定を省略）を指定してもできるけど．．．。
cpt_nile <- cpt.meanvar(Nile, method='AMOC')
# cpt_nile <- cpt.meanvar(Nile)
plot(cpt_nile)
参考までにパラメータのいじり方も少し練習しておこうか。
ラジャー。
変化点の数は、ペナルティの設定を変更すると変化するよ。具体的には引数penaltyにManualを指定して、引数pen.valueの値を変えることで変化点の数が変わるよ。pen.valueは小さな値にすると変化点が増えて、大きな値にすると変化点が減るよ。
par(mfrow=c(3, 2))
cpt_nile_pv1 <- cpt.meanvar(Nile, method='PELT', penalty='Manual', pen.value=1)
plot(cpt_nile_pv1, main='pen.value=1')
cpt_nile_pv5 <- cpt.meanvar(Nile, method='PELT', penalty='Manual', pen.value=5)
plot(cpt_nile_pv5, main='pen.value=5')
cpt_nile_pv10 <- cpt.meanvar(Nile, method='PELT', penalty='Manual', pen.value=10)
plot(cpt_nile_pv10, main='pen.value=10')
cpt_nile_pv20 <- cpt.meanvar(Nile, method='PELT', penalty='Manual', pen.value=20)
plot(cpt_nile_pv20, main='pen.value=20')
cpt_nile_pv50 <- cpt.meanvar(Nile, method='PELT', penalty='Manual', pen.value=50)
plot(cpt_nile_pv50, main='pen.value=50')
cpt_nile_pv100 <- cpt.meanvar(Nile, method='PELT', penalty='Manual', pen.value=100)
plot(cpt_nile_pv100, main='pen.value=100')
par(mfrow=c(1, 1))
なるほど。でも、よさげなpen.valの値はどうやって見つけるの？
完璧な方法は無さそうだけど、自動探索機能と探索したプロット結果を見れば、ある程度妥当なpen.val値を求めることができるみたい。
pen.valを自動探索するには、penaltyにCROPSを指定、pen.valはベクトルで検索範囲（上下限値）を指定すればOKだよ。
# 妥当な変化点数、ペナルティを求めるための自動探索
cpt_nile_search <- cpt.meanvar(Nile, method='PELT', penalty='CROPS', pen.value=c(10, 500))
# 検出された変化点
print(cpts.full(cpt_nile_search))
# 対応するペナルティ（pen.valの値）
print(pen.value.full(cpt_nile_search))
自動探索の診断結果をプロットするには、plot()関数の引数にdiagnostic=TRUEを指定してね。
# 診断プロット
plot(cpt_nile_search, diagnostic=TRUE)
最適な変化の点の数は、ペナルティの変化量が小さくなるところ。曲線データの肘に相当するあたりがよさげらしいよ。
肘？がどこだか明確には分からないけど、なんとなく変化点数が4か3あたりかな？
探索結果の情報をもとにグラフを描くだけだったら、引数ncptsに変化点の数を指定してプロットすればOKだよ。ただし探索をしていない変化点の数値を指定するとエラーになっちゃうから注意ね（今回の場合は11,8,4,3,1から選択できる）。
par(mfrow=c(3, 2))
plot(cpt_nile_search, ncpts=11, main='ncpts=11')
plot(cpt_nile_search, ncpts=8, main='ncpts=8')
plot(cpt_nile_search, ncpts=4, main='ncpts=4')
plot(cpt_nile_search, ncpts=3, main='ncpts=3')
plot(cpt_nile_search, ncpts=1, main='ncpts=1')
plot(cpt_nile_search, diagnostic=TRUE)
par(mfrow=c(1, 1))
探索結果から求められたpen.valueの値をつかって、再計算を実行してもいいよ。たとえば、探索で見つかった変化点4個を再現したい場合はこんな感じで。
cpts.full(cpt_nile_search)[3,] # 探索結果の3番目が変化点を4個検出
pv <- pen.value.full(cpt_nile_search)[3] # 探索結果3番目のデータのpen.value値
pv
cpt_nile_cp4 <- cpt.meanvar(Nile, method='PELT', penalty='Manual', pen.value=pv)
plot(cpt_nile_cp4, main='pen.value=11.73')
試しに．．．
Sepal_Length <- iris[, 1] # iris Sepal Length
Sepal_Width <- iris[, 2] # iris Sepal Width
Petal_Length <- iris[, 3] # iris Petal Length
Petal_Width <- iris[, 4] # iris Petal Width
par(mfrow=c(2, 2))
cpt_Sepal_Length <- cpt.meanvar(Sepal_Length, method='PELT', penalty='Manual', pen.value=40)
plot(cpt_Sepal_Length, cpt.width=1, xlab='iris index', ylab='value')
title('Iris Sepal Length')
cpt_Sepal_Width <- cpt.meanvar(Sepal_Width, method='PELT', penalty='Manual', pen.value=50)
plot(cpt_Sepal_Width, cpt.width=1, xlab='iris index', ylab='value')
title('Iris Sepal Width')
cpt_Petal_Length <- cpt.var(Petal_Length, method='PELT')
plot(cpt_Petal_Length, cpt.width=1, xlab='iris index', ylab='value')
title('Iris Petal Length')
cpt_Petal_Width <- cpt.var(Petal_Width, method='PELT')
plot(cpt_Petal_Width, cpt.width=1, xlab='iris index', ylab='value')
title('Iris Petal Width')
par(mfrow=c(1, 1))
Irisデータでchangepointを使ってみた。3種類の花のデータが順番（setosa,versicolor,virginica）に並んでいるから、花の種類が変わるところで変化点を検出できるかもと．．．。
setosa(1-50)は特徴が大きく違うから検出も簡単。versicolor(51-100)とvirginica(101-150)はよく似ているから、特徴量Sepalでは変化点がうまく検出できなかったみたい。
無理くりIris？使い方としてはちょっとどうかなぁ．．．。まぁ練習なのでヨシとしよう。
今回のchangepointを使った変化点検出は以上で。またね！
お疲れ様でした。
0 件のコメント :
コメントを投稿