読者です 読者をやめる 読者になる 読者になる

でたぁっ 感動と失敗の備忘録

データ解析を担当することになったサラリーマンの備忘録

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

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
f:id:deta:20130514053239p:plain

change in mean
f:id:deta:20130514053302p:plain

change in mean and variance
f:id:deta:20130514053316p:plain

といった感じです。具体的なデータではどのようになるの?ということで日経平均株価を試してみましたのでその結果を紹介します。データはダウンロードサイトから 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 

f:id:deta:20130514054335p:plain

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


参考
changepoint: An R Package for Changepoint Analysis
データマイニングによる異常検知
データマイニングによる異常検知
異常検知(変化点検出)のパッケージを作ってみた
隠れマルコフモデルで異常検知(R Advent Calendar2012)