R changepoint パッケージちょっと試した

　よく見ている R-bloggers で R で changepoint パッケージを使ったコードが掲載されていた。（実際のサイトはこちら）

変化点検出は、これまでもいくつか書籍やブログ、R勉強会のTokyoRなどで紹介されている。

・データマイニングによる異常検知　ChangeFinderが掲載されている山西さんの書籍
・異常検知（変化点検出）のパッケージを作ってみた TokyoR主催者yokkunsさんの日記
・隠れマルコフモデルで異常検知（R Advent Calendar2012） teramonagiさんのブログ（コードはこちら）

changepoint パッケージは初めてだったので、ちょっとだけ試してみた。まずは inside-R にあったサンプルコードを実行。サンプルデータは正規分布の乱数です。

library(changepoint)

# change in variance
x=c(rnorm(100,0,1),rnorm(100,0,10))
ansvar=cpt.var(x)
plot(ansvar)
print(ansvar)
 
# change in mean
y=c(rnorm(100,0,1),rnorm(100,5,1))
ansmean=cpt.mean(y)
plot(ansmean,cpt.col='blue')
print(ansmean)
 
# change in mean and variance
z=c(rnorm(100,0,1),rnorm(100,2,10))
ansmeanvar=cpt.meanvar(z)
plot(ansmeanvar,cpt.width=3)
print(ansmeanvar)

マニュアルにも cpt.mean、cpt.var、cpt.meanvar を見てから始めて．．．といった説明もあるので、基本的な関数はこの3つ？　それぞれの plot はこのようになります。

change in variance

change in mean

change in mean and variance

といった感じです。具体的なデータではどのようになるの？ということで日経平均株価を試してみましたのでその結果を紹介します。データはダウンロードサイトから csv ファイルを準備しました。（R では yokkunsさんが開発した株価が取得できる RFinanceYJ というパッケージもあります）csv ファイルは5列目が終値です。また、関数は cpt.meanvar、引数の method は AMOC、PELT、SegNeigh、BinSeg の中から PELT を選択しております。

library(changepoint)

## ---日経平均日足データの読込み
setwd("d:/R/changepoint") # 作業ディレクトリの変更

d <- read.csv("日経平均日足.csv", header=T) # データフレームで読込み

head(d)
       日付    始値    高値    安値    終値
1  2012/5/8 9190.50 9207.56 9159.47 9181.65
2  2012/5/9 9112.72 9115.94 9021.20 9045.06
3 2012/5/10 9013.26 9075.63 8985.90 9009.65
4 2012/5/11 9019.40 9050.61 8944.63 8953.31
5 2012/5/14 8986.22 9031.09 8947.82 8973.84
6 2012/5/15 8910.85 8930.78 8838.78 8900.74


ansmeanvar=cpt.meanvar(d[,5], method='PELT')
plot(ansmeanvar,cpt.width=3)
print(ansmeanvar)

Class 'cpt' : Changepoint Object
       ~~   : S4 class containing 10 slots with names
              date data.set cpttype method test.stat pen.type pen.value cpts ncpts.max param.est 

Created on  : Sun May 12 17:57:17 2013 

summary(.)  :
----------
Changepoint type      : Change in mean and variance 
Method of analysis    : PELT 
Test Statistic  : Normal 
Type of penalty       : SIC with value, 11.04292 
Maximum no. of cpts   : Inf 
Number of changepoints: 42 

一応変化点が検出されてっるぽい。まだまだ勉強不足なので時間があったらもう少し調べる。

参考
・changepoint: An R Package for Changepoint Analysis
・データマイニングによる異常検知

・異常検知（変化点検出）のパッケージを作ってみた
・隠れマルコフモデルで異常検知（R Advent Calendar2012）