#title Change Point
[[TableOfContents]]

slope가 변하는 지점을 찾는다.

==== 변경점 ====
{{{
library(segmented)
options(warn = -1)

change_point <- function(y, chart=F, n=10){
  #y <- sin(sort(runif(100,0,5)))*10+100
  #n=5
  #chart=T
  
  x <- 1:length(y)
  #x <- (1:length(y)) * (max(y) - min(y)) / length(y)

  for(cp_cnt in 1:n){
    fit_lm = lm(y ~ x)  # intercept-only model
    fit_tmp = try(segmented(fit_lm, seg.Z = ~x, npsi = cp_cnt),silent=T)
    
    if (class(fit_tmp)[1] == "try-error"){
      break
      cp_cnt <- cp_cnt - 1
    } else if (length(fit_tmp$psi) == 0){
      break
      cp_cnt <- cp_cnt - 1
    }
    
    r2 <- cor(y, fitted(fit_tmp))^2  
    if(r2 >= 0.9){
      fit_segmented <- fit_tmp
      break
    }
    
    fit_segmented <- fit_tmp
  }
  
  if(1==chart){
    plot(fit_segmented)
    points(x,y)
    lines.segmented(fit_segmented)
    points.segmented(fit_segmented)  
  }
  
  cp <- fit_segmented$psi[,2]
  return(cp)
}


set.seed(1024)
y <- sin(sort(runif(100,0,5)))*10+100 + sin(sort(runif(50,0,5)))*10+100
cp <- change_point(y, T)
}}}

결과
attachment:ChangePoint/change_point_2.png


==== 최근 변경점 ====
{{{
library(segmented)
options(warn = -1)

recent_change_point <- function(y, chart=F, n=10){
  #y <- sin(sort(runif(100,0,5)))*10+100
  #n=5
  #chart=T
  
  x <- 1:length(y)
  #x <- (1:length(y)) * (max(y) - min(y)) / length(y)

  for(cp_cnt in 1:n){
    fit_lm = lm(y ~ x)  # intercept-only model
    fit_tmp = try(segmented(fit_lm, seg.Z = ~x, npsi = cp_cnt),silent=T)
    
    if (class(fit_tmp) == "try-error"){
      break
      cp_cnt <- cp_cnt - 1
    } else if (length(fit_tmp$psi) == 0){
      break
      cp_cnt <- cp_cnt - 1
    }
    
    r2 <- cor(y, fitted(fit_tmp))^2  
    if(r2 >= 0.9){
      fit_segmented <- fit_tmp
      break
    }
    
    fit_segmented <- fit_tmp
  }
  
  if(1==chart){
    plot(fit_segmented)
    points(x,y)
    lines.segmented(fit_segmented)
    points.segmented(fit_segmented)  
  }
  
  cp <- fit_segmented$psi[,2]
  
  start_x <- round(tail(cp,1), 2)
  y2 <- y[x >= start_x]
  x2 <- x[x >= start_x]
  return(data.frame(x=x2, y=y2))
}


set.seed(1024)
y <- sin(sort(runif(100,0,5)))*10+100 + sin(sort(runif(50,0,5)))*10+100
x <- 1:length(y)
plot(x, y)
lines(recent_change_point(y))
}}}

결과
attachment:ChangePoint/recent_change_point_3.png