관리자.. 2012. 11. 28. 06:26
오늘은 bfast 예제 실행 해 보고, 내 데이터에 적용하고 그러다 낮이 갔다. 수업 들어가서 조교노릇 하고, 여튼 괜찮았고 생산적이었음. 
R에서 시계열 자료의 구조적 변화를 찾는 패키지가 몇 개 있는데, 그 중에 우리 분야에서 자주 언급되는 bfast 랑 bise, 그리고 새로 소개된 몇 개 season 이라던가 보다 범용 패키지도 시험 해 보고 했다.
bfast 에서는 break point를 감지하는데, 이게 그리 정확하게 되진 않았다. 당연히 break point 의 정의에 관련되는 파라미터를 어떻게 조정하느냐에 결과가 많이 영향을 받는데, 그건 괜찮다 치고, 지난 번에 시계열 수업에서 제기됐던 문제가 여기도 있었다. 경향이 강할 때 계절성을 찾는 것, 혹은 계절성이 강할 때 경향을 찾는 것이 까다롭다는 것. 그래도 예전에 쓰던, 오히려 상당히 진보된 방법론인 SSA 대신 이걸 써 보는 이유는, SSA가 국지적으로 최적화 하면서 눈에 보이는 패턴을 잘 찾아내긴 하지만 외려 전역적인 특성을 기술하는 데에는 지나치게 복잡하다는 것, 그래서 어느 선에서 그어야 한다면 차라리 고전적인 계절성 추출이 더 내 연구 목적에 맞지 않냐, 그런 문제.

패키지의 기본 예제를 좀 보자면,  
## Simulated Data
plot(simts) # stl object containing simulated NDVI time series

datats <- ts(rowSums(simts$time.series)) # sum of all the components (season,abrupt,remainder)
tsp(datats) <- tsp(simts$time.series) # assign correct time series attributes
plot(datats)

fit <- bfast(datats,h=0.15, season="dummy", max.iter=1) 
plot(fit,sim=simts)
fit # prints out whether breakpoints are detected in the seasonal and trend component

## Real data 이 부분을 보자
## The data should be a regular ts() object without NA's. R 기본 ts 객체면 됨.
## See Fig. 8 b in reference
plot(harvest, ylab="NDVI") # MODIS 16-day cleaned and interpolated NDVI time series 


16일 간격 NDVI 자료. 2000년 부터 2009년 까지 이어놓으니 변화가 한 눈에 보인다. 이 정도 변화가 보인다는 건 토지피복이 2995년 전후해서 크게 변했다는 것.


(rdist <- 10/length(harvest)) # ratio of distance between breaks (time steps) and length of the time series fit <- bfast(harvest,h=rdist, season="harmonic", max.iter=1,breaks=2)
이 부분에서 break 가 2 개 일 거라고 지정을 해 줘야 했다 이 방법론에서는. 이걸 바꿔가면서 테스트 해 보는 게 의미있었다.
plot(fit)

fit 된 결과를 도시 해 보면, 2 개의 급격한 변화 지점을 찾으라는 명령 대로 3개의 구간으로 나눠 놓은 것을 볼 수 있다. 세 번째 그림이 트렌드고,


## plot anova and slope of the trend identified trend segments plot(fit, ANOVA=TRUE) ## plot the trend component and identify the break with the largest magnitude of change plot(fit,type="trend",largest=TRUE) ## plot all the different available plots 
plot(fit,type="all") 


작성 중 
졸려서 집에 이란 가야겠다;;
--
이거 말고도 작성하다 엎어진 R 일기가 몇 개 있는데, 여튼 올 해 가기 전에 다 정리하겠다는 의지로 @_@
bfast 는 앞으로 이삼일 계속 해야니까 금방 완성할 수 있겠고, 아래에 보면 PLFA 분석도 오늘 마르코 와서 또 한 번 했는데, 마르코씨 컴퓨터에서 작업을 하니까 아무래도 정리가 잘 안된다 마치고 난 후에. 완전히 끝나면 꼭꼭 정리하겠음 (주성분 분석 좀 사례 들어서 세세히 보고자 함). BCART는 사실 내가 좀 보다 밀어놔서 시간이 더 걸릴 듯..  연말 휴가 땐 R 해야지 히힛.